[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]

details vigra/matrix.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*        Copyright 2004 by Gunnar Kedenburg and Ullrich Koethe         */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    ( Version 1.3.2, Jan 27 2005 )                                    */
00008 /*    You may use, modify, and distribute this software according       */
00009 /*    to the terms stated in the LICENSE file included in               */
00010 /*    the VIGRA distribution.                                           */
00011 /*                                                                      */
00012 /*    The VIGRA Website is                                              */
00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00014 /*    Please direct questions, bug reports, and contributions to        */
00015 /*        koethe@informatik.uni-hamburg.de                              */
00016 /*                                                                      */
00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
00020 /*                                                                      */
00021 /************************************************************************/
00022 
00023 
00024 #ifndef VIGRA_MATRIX_HXX
00025 #define VIGRA_MATRIX_HXX
00026 
00027 #include <cmath>
00028 #include <iosfwd>
00029 #include <iomanip>
00030 #include "vigra/multi_array.hxx"
00031 #include "vigra/mathutil.hxx"
00032 
00033 
00034 namespace vigra
00035 {
00036 
00037 namespace linalg
00038 {
00039 
00040 template <class T, class C>
00041 inline std::size_t rowCount(const MultiArrayView<2, T, C> &x);
00042 
00043 template <class T, class C>
00044 inline std::size_t columnCount(const MultiArrayView<2, T, C> &x);
00045 
00046 template <class T, class C>
00047 MultiArrayView <2, T, C>
00048 rowVector(MultiArrayView <2, T, C> const & m, int d);
00049 
00050 template <class T, class C>
00051 MultiArrayView <2, T, C>
00052 columnVector(MultiArrayView<2, T, C> const & m, int d);
00053 
00054 template <class T, class C>
00055 T squaredNorm(const MultiArrayView<2, T, C> &a);
00056 
00057 template <class T, class ALLOC>
00058 class TemporaryMatrix;
00059 
00060 template <class T, class C1, class C2>
00061 void transpose(const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &r);
00062 
00063 template <class T, class C>
00064 bool isSymmetric(const MultiArrayView<2, T, C> &v);
00065 
00066 
00067 /********************************************************/
00068 /*                                                      */
00069 /*                        Matrix                        */
00070 /*                                                      */
00071 /********************************************************/
00072 
00073 /** Matrix class.
00074 
00075     This is the basic class for all linear algebra computations. Matrices are
00076     strored in a <i>column-major</i> format, i.e. the row index is varying fastest.
00077     This is the same format as in the lapack and gmm++ libraries, so it will
00078     be easy to interface these libraries. In fact, if you need optimized
00079     high performance code, you should use them. The VIGRA linear algebra
00080     functionality is provided for smaller problems and rapid prototyping
00081     (no one wants to spend half the day installing a new library just to 
00082     discover that the new algorithm idea didn't work anyway).
00083 
00084     <b>See also:</b>
00085     <ul>
00086     <li> \ref LinearAlgebraFunctions
00087     </ul>
00088 
00089     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00090     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00091         Namespaces: vigra and vigra::linalg
00092 */
00093 template <class T, class ALLOC = std::allocator<T> >
00094 class Matrix
00095 : public MultiArray<2, T, ALLOC>
00096 {
00097     typedef MultiArray<2, T, ALLOC> BaseType;
00098     
00099   public:
00100     typedef Matrix<T, ALLOC>                        matrix_type;
00101     typedef TemporaryMatrix<T, ALLOC>               temp_type;
00102     typedef MultiArrayView<2, T, UnstridedArrayTag> view_type;
00103     typedef typename BaseType::value_type           value_type;
00104     typedef typename BaseType::pointer              pointer;
00105     typedef typename BaseType::const_pointer        const_pointer;
00106     typedef typename BaseType::reference            reference;
00107     typedef typename BaseType::const_reference      const_reference;
00108     typedef typename BaseType::difference_type      difference_type;
00109     typedef ALLOC                                   allocator_type;
00110     
00111         /** default constructor
00112          */
00113     Matrix() 
00114     {}
00115 
00116         /** construct with given allocator
00117          */
00118     explicit Matrix(ALLOC const & alloc)
00119     : BaseType(alloc)
00120     {}
00121 
00122         /** construct with given shape and init with all 
00123             elements with zero. Note that the order of the axes is
00124             <tt>difference_type(rows, columns)</tt> which
00125             is the opposite of the usual VIGRA convention.
00126          */
00127     explicit Matrix(const difference_type &shape,
00128                     ALLOC const & alloc = allocator_type())
00129     : BaseType(shape, alloc)
00130     {}
00131 
00132         /** construct with given shape and init with all 
00133             elements with zero. Note that the order of the axes is
00134             <tt>(rows, columns)</tt> which
00135             is the opposite of the usual VIGRA convention.
00136          */
00137     Matrix(std::size_t rows, std::size_t columns,
00138                     ALLOC const & alloc = allocator_type())
00139     : BaseType(difference_type(rows, columns), alloc)
00140     {}
00141 
00142         /** construct with given shape and init with all 
00143             elements with the constant \a init. Note that the order of the axes is
00144             <tt>difference_type(rows, columns)</tt> which
00145             is the opposite of the usual VIGRA convention.
00146          */
00147     Matrix(const difference_type &shape, const_reference init,
00148            allocator_type const & alloc = allocator_type())
00149     : BaseType(shape, init, alloc)
00150     {}
00151 
00152         /** construct with given shape and init with all 
00153             elements with the constant \a init. Note that the order of the axes is
00154             <tt>(rows, columns)</tt> which
00155             is the opposite of the usual VIGRA convention.
00156          */
00157     Matrix(std::size_t rows, std::size_t columns, const_reference init,
00158            allocator_type const & alloc = allocator_type())
00159     : BaseType(difference_type(rows, columns), init, alloc)
00160     {}
00161 
00162         /** construct with given shape and copy data from C-style array \a init.
00163             Data in this array are expected to be given in column-major
00164             order (the C standard order) and will automatically be
00165             converted to the required column-major format. Note that the order of the axes is
00166             <tt>difference_type(rows, columns)</tt> which
00167             is the opposite of the usual VIGRA convention.
00168          */
00169     Matrix(const difference_type &shape, const_pointer init,
00170            allocator_type const & alloc = allocator_type())
00171     : BaseType(shape, alloc)
00172     {
00173         difference_type trans(shape[1], shape[0]);
00174         linalg::transpose(MultiArrayView<2, T>(trans, const_cast<pointer>(init)), *this);
00175     }
00176 
00177         /** construct with given shape and copy data from C-style array \a init.
00178             Data in this array are expected to be given in column-major
00179             order (the C standard order) and will automatically be
00180             converted to the required column-major format. Note that the order of 
00181             the axes is <tt>(rows, columns)</tt> which
00182             is the opposite of the usual VIGRA convention.
00183          */
00184     Matrix(std::size_t rows, std::size_t columns, const_pointer init,
00185            allocator_type const & alloc = allocator_type())
00186     : BaseType(difference_type(rows, columns), alloc)
00187     {
00188         difference_type trans(columns, rows);
00189         linalg::transpose(MultiArrayView<2, T>(trans, const_cast<pointer>(init)), *this);
00190     }
00191 
00192         /** copy constructor. Allocates new memory and 
00193             copies tha data.
00194          */
00195     Matrix(const Matrix &rhs)
00196     : BaseType(rhs)
00197     {}
00198 
00199         /** construct from temporary matrix, which looses its data.
00200             
00201             This operation is equivalent to
00202             \code
00203                 TemporaryMatrix<T> temp = ...;
00204                 
00205                 Matrix<T> m;
00206                 m.swap(temp);
00207             \endcode
00208          */
00209     Matrix(const TemporaryMatrix<T, ALLOC> &rhs)
00210     : BaseType(rhs.allocator())
00211     {
00212         this->swap(const_cast<TemporaryMatrix<T, ALLOC> &>(rhs));
00213     }
00214     
00215         /** construct from a MultiArrayView. Allocates new memory and 
00216             copies tha data. \a rhs is assumed to be in column-major order already.
00217          */
00218     template<class U, class C>
00219     Matrix(const MultiArrayView<2, U, C> &rhs)
00220     : BaseType(rhs)
00221     {}
00222 
00223         /** assignment.
00224             If the size of \a rhs is the same as the matrix's old size, only the data
00225             are copied. Otherwise, new storage is allocated, which invalidates 
00226             all objects (array views, iterators) depending on the matrix.
00227          */
00228     Matrix & operator=(const Matrix &rhs)
00229     {
00230         BaseType::operator=(rhs); // has the correct semantics already
00231         return *this;
00232     }
00233 
00234         /** assign a temporary matrix. This is implemented by swapping the data
00235             between the two matrices, so that all depending objects 
00236             (array views, iterators) ar invalidated.
00237          */
00238     Matrix & operator=(const TemporaryMatrix<T, ALLOC> &rhs)
00239     {
00240         this->swap(const_cast<TemporaryMatrix<T, ALLOC> &>(rhs));
00241         return *this;
00242     }
00243 
00244         /** assignment from arbitrary 2-dimensional MultiArrayView.<br>
00245             If the size of \a rhs is the same as the matrix's old size, only the data
00246             are copied. Otherwise, new storage is allocated, which invalidates 
00247             all objects (array views, iterators) depending on the matrix. 
00248             \a rhs is assumed to be in column-major order already.
00249          */
00250     template <class U, class C>
00251     Matrix & operator=(const MultiArrayView<2, U, C> &rhs)
00252     {
00253         BaseType::operator=(rhs); // has the correct semantics already
00254         return *this;
00255     }
00256     
00257         /** Create a matrix view that represents the row vector of row \a d.
00258          */
00259     view_type rowVector(std::size_t d) const
00260     {
00261         return vigra::linalg::rowVector(*this, d);
00262     }
00263     
00264         /** Create a matrix view that represents the column vector of column \a d.
00265          */
00266     view_type columnVector(std::size_t d) const
00267     {
00268         return vigra::linalg::columnVector(*this, d);
00269     }
00270     
00271         /** number of rows (height) of the matrix.
00272         */
00273     std::size_t rowCount() const
00274     {
00275         return this->m_shape[0];
00276     }
00277     
00278         /** number of columns (width) of the matrix.
00279         */
00280     std::size_t columnCount() const
00281     {
00282         return this->m_shape[1];
00283     }
00284     
00285         /** number of elements (width*height) of the matrix.
00286         */
00287     std::size_t elementCount() const
00288     {
00289         return rowCount()*columnCount();
00290     }
00291        
00292         /** check whether the matrix is symmetric.
00293         */
00294     bool isSymmetric() const
00295     {
00296         return vigra::linalg::isSymmetric(*this);
00297     }
00298     
00299 #ifdef DOXYGEN 
00300 // repeat the index functions for documentation. In real code, they are inherited.
00301 
00302         /** read/write access to matrix element <tt>(row, column)</tt>.
00303             Note that the order of the argument is the opposite of the usual 
00304             VIGRA convention due to column-major matrix order.
00305         */
00306     value_type & operator()(std::size_t row, std::size_t column);
00307 
00308         /** read access to matrix element <tt>(row, column)</tt>.
00309             Note that the order of the argument is the opposite of the usual 
00310             VIGRA convention due to column-major matrix order.
00311         */
00312     value_type operator()(std::size_t row, std::size_t column) const;
00313 #endif
00314 
00315         /** squared Frobenius norm. Sum of squares of the matrix elements.
00316         */
00317     value_type squaredNorm() const
00318     {
00319         return vigra::linalg::squaredNorm(*this);
00320     }
00321     
00322         /** Frobenius norm. Root of sum of squares of the matrix elements.
00323         */
00324     value_type norm() const
00325     {
00326         return VIGRA_CSTD::sqrt(squaredNorm());
00327     }
00328     
00329         /** transpose matrix in-place (precondition: matrix must be square)
00330          */
00331     Matrix & transpose();
00332     
00333         /** add \a other to this (sizes must match).
00334          */
00335     template <class U, class C>
00336     Matrix & operator+=(MultiArrayView<2, U, C> const & other);
00337     
00338         /** subtract \a other from this (sizes must match).
00339          */
00340     template <class U, class C>
00341     Matrix & operator-=(MultiArrayView<2, U, C> const & other);
00342     
00343         /** scalar multiply this with \a other
00344          */
00345     Matrix & operator*=(T other);
00346     
00347         /** scalar devide this by \a other
00348          */
00349     Matrix & operator/=(T other);
00350 };
00351 
00352 template <class T, class ALLOC>
00353 Matrix<T, ALLOC> & Matrix<T, ALLOC>::transpose()
00354 {
00355     const unsigned int cols = columnCount();
00356     vigra_precondition(cols == rowCount(),
00357         "Matrix::transpose(): in-place transposition requires square matrix.");
00358     for(unsigned int i = 0; i < cols; ++i)
00359         for(unsigned int j = i+1; j < cols; ++j)
00360             std::swap((*this)(j, i), (*this)(i, j));
00361     return *this;
00362 }
00363 
00364 template <class T, class ALLOC>
00365 template <class U, class C>
00366 Matrix<T, ALLOC> & Matrix<T, ALLOC>::operator+=(MultiArrayView<2, U, C> const & other)
00367 {
00368     const unsigned int rows = rowCount();
00369     const unsigned int cols = columnCount();
00370     vigra_precondition(rows == vigra::linalg::rowCount(other) && cols == vigra::linalg::columnCount(other),
00371        "Matrix::operator+=(): Shape mismatch.");
00372     
00373     for(unsigned int i = 0; i < cols; ++i)
00374         for(unsigned int j = 0; j < rows; ++j)
00375             (*this)(j, i) += other(j, i);
00376     return *this;
00377 }
00378 
00379 template <class T, class ALLOC>
00380 template <class U, class C>
00381 Matrix<T, ALLOC> & Matrix<T, ALLOC>::operator-=(MultiArrayView<2, U, C> const & other)
00382 {
00383     const unsigned int rows = rowCount();
00384     const unsigned int cols = columnCount();
00385     vigra_precondition(rows == vigra::linalg::rowCount(other) && cols == vigra::linalg::columnCount(other),
00386        "Matrix::operator-=(): Shape mismatch.");
00387     
00388     for(unsigned int i = 0; i < cols; ++i)
00389         for(unsigned int j = 0; j < rows; ++j)
00390             (*this)(j, i) -= other(j, i);
00391     return *this;
00392 }
00393 
00394 template <class T, class ALLOC>
00395 Matrix<T, ALLOC> & Matrix<T, ALLOC>::operator*=(T other)
00396 {
00397     const unsigned int rows = rowCount();
00398     const unsigned int cols = columnCount();
00399     
00400     for(unsigned int i = 0; i < cols; ++i)
00401         for(unsigned int j = 0; j < rows; ++j)
00402             (*this)(j, i) *= other;
00403     return *this;
00404 }
00405 
00406 template <class T, class ALLOC>
00407 Matrix<T, ALLOC> & Matrix<T, ALLOC>::operator/=(T other)
00408 {
00409     const unsigned int rows = rowCount();
00410     const unsigned int cols = columnCount();
00411     
00412     for(unsigned int i = 0; i < cols; ++i)
00413         for(unsigned int j = 0; j < rows; ++j)
00414             (*this)(j, i) /= other;
00415     return *this;
00416 }
00417 
00418 // TemporaryMatrix is provided as an optimization: Functions returning a matrix can 
00419 // use TemporaryMatrix to make explicit that it was allocated as a temporary data structure.
00420 // Functions receiving a TemporaryMatrix can thus often avoid to allocate new temporary 
00421 // memory.
00422 template <class T, class ALLOC = std::allocator<T> >
00423 class TemporaryMatrix
00424 : public Matrix<T, ALLOC>
00425 {
00426     typedef Matrix<T, ALLOC> BaseType;
00427   public:
00428     typedef Matrix<T, ALLOC>                        matrix_type;
00429     typedef TemporaryMatrix<T, ALLOC>               temp_type;
00430     typedef MultiArrayView<2, T, UnstridedArrayTag> view_type;
00431     typedef typename BaseType::value_type           value_type;
00432     typedef typename BaseType::pointer              pointer;
00433     typedef typename BaseType::const_pointer        const_pointer;
00434     typedef typename BaseType::reference            reference;
00435     typedef typename BaseType::const_reference      const_reference;
00436     typedef typename BaseType::difference_type      difference_type;
00437     typedef ALLOC                                   allocator_type;
00438 
00439     TemporaryMatrix(std::size_t rows, std::size_t columns)
00440     : BaseType(rows, columns, ALLOC())
00441     {}
00442 
00443     TemporaryMatrix(std::size_t rows, std::size_t columns, const_reference init)
00444     : BaseType(rows, columns, init, ALLOC())
00445     {}
00446 
00447     template<class U, class C>
00448     TemporaryMatrix(const MultiArrayView<2, U, C> &rhs)
00449     : BaseType(rhs)
00450     {}
00451     
00452     TemporaryMatrix(const TemporaryMatrix &rhs)
00453     : BaseType()
00454     {
00455         this->swap(const_cast<TemporaryMatrix &>(rhs));
00456     }
00457     
00458     TemporaryMatrix & transpose()
00459     {
00460         BaseType::transpose();
00461         return *this;
00462     }
00463     
00464     template <class U, class C>
00465     TemporaryMatrix & operator+=(MultiArrayView<2, U, C> const & other)
00466     {
00467         BaseType::operator+=(other);
00468         return *this;
00469     }
00470     
00471     template <class U, class C>
00472     TemporaryMatrix & operator-=(MultiArrayView<2, U, C> const & other)
00473     {
00474         BaseType::operator-=(other);
00475         return *this;
00476     }
00477 
00478     TemporaryMatrix & operator*=(T other)
00479     {
00480         BaseType::operator*=(other);
00481         return *this;
00482     }
00483     
00484     TemporaryMatrix & operator/=(T other)
00485     {
00486         BaseType::operator/=(other);
00487         return *this;
00488     }
00489   private:
00490 
00491     TemporaryMatrix &operator=(const TemporaryMatrix &rhs); // not implemented
00492 };
00493 
00494 
00495 /** \addtogroup LinearAlgebraFunctions Matrix functions
00496  */
00497 //@{
00498  
00499     /** Number of rows of a matrix represented as a <tt>MultiArrayView&lt;2,...&gt;</tt>
00500     
00501     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00502     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00503         Namespaces: vigra and vigra::linalg
00504      */ 
00505 template <class T, class C>
00506 inline std::size_t rowCount(const MultiArrayView<2, T, C> &x)
00507 {
00508     return x.shape(0);
00509 }
00510 
00511     /** Number of columns of a matrix represented as a <tt>MultiArrayView&lt;2,...&gt;</tt>
00512     
00513     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00514     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00515         Namespaces: vigra and vigra::linalg
00516      */ 
00517 template <class T, class C>
00518 inline std::size_t columnCount(const MultiArrayView<2, T, C> &x)
00519 {
00520     return x.shape(1);
00521 }
00522 
00523     /** Create a row vector view for row \a d of the matrix \a m
00524     
00525     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00526     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00527         Namespaces: vigra and vigra::linalg
00528      */ 
00529 template <class T, class C>
00530 MultiArrayView <2, T, C>
00531 rowVector(MultiArrayView <2, T, C> const & m, int d)
00532 {
00533     typedef typename MultiArrayView <2, T, C>::difference_type Shape;
00534     return m.subarray(Shape(d, 0), Shape(d+1, columnCount(m)));
00535 }
00536 
00537     /** Create a column vector view for column \a d of the matrix \a m
00538     
00539     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00540     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00541         Namespaces: vigra and vigra::linalg
00542      */ 
00543 template <class T, class C>
00544 MultiArrayView <2, T, C>
00545 columnVector(MultiArrayView<2, T, C> const & m, int d)
00546 {
00547     typedef typename MultiArrayView <2, T, C>::difference_type Shape;
00548     return m.subarray(Shape(0, d), Shape(rowCount(m), d+1));
00549 }
00550 
00551     /** Check whether matrix \a m is symmetric.
00552     
00553     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00554     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00555         Namespaces: vigra and vigra::linalg
00556      */ 
00557 template <class T, class C>
00558 bool
00559 isSymmetric(MultiArrayView<2, T, C> const & m)
00560 {
00561     const unsigned int size = rowCount(m);
00562     if(size != columnCount(m))
00563         return false;
00564         
00565     for(unsigned int i = 0; i < size; ++i)
00566         for(unsigned int j = i+1; j < size; ++j)
00567             if(m(j, i) != m(i, j))
00568                 return false;
00569     return true;
00570 }
00571 
00572 
00573     /** calculate the squared Frobenius norm of a matrix. 
00574         Equal to the sum of squares of the matrix elements.
00575     
00576     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00577     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00578         Namespaces: vigra and vigra::linalg
00579      */ 
00580 template <class T, class C>
00581 T squaredNorm(const MultiArrayView<2, T, C> &a)
00582 {
00583     const unsigned int rows = rowCount(a);
00584     const unsigned int cols = columnCount(a);
00585     T ret = NumericTraits<T>::zero();
00586     for(unsigned int j = 0; j < cols; ++j)
00587         for(unsigned int i = 0; i < rows; ++i)
00588             ret += sq(a(i, j));
00589     return ret;
00590 }
00591 
00592     /** calculate the squared norm of a vector. 
00593         Equal to the sum of squares of the vector elements.
00594     
00595     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00596     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00597         Namespaces: vigra and vigra::linalg
00598      */ 
00599 template <class T, class C>
00600 T squaredNorm(const MultiArrayView<1, T, C> &a)
00601 {
00602     const unsigned int size = a.elementCount();
00603     T ret = NumericTraits<T>::zero();
00604     for(unsigned int i = 0; i < size; ++i)
00605         ret += sq(a(i));
00606     return ret;
00607 }
00608 
00609     /** calculate the Frobenius norm of a matrix or vector. 
00610         Equal to the square root of sum of squares of the matrix elements.
00611     
00612     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00613     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00614         Namespaces: vigra and vigra::linalg
00615      */ 
00616 template <unsigned int N, class T, class C>
00617 T norm(const MultiArrayView<N, T, C> &a)
00618 {
00619     return VIGRA_CSTD::sqrt(squaredNorm(a));
00620 }
00621 
00622     /** initialize the given square matrix as an identity matrix.
00623     
00624     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00625     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00626         Namespaces: vigra and vigra::linalg
00627      */ 
00628 template <class T, class C>
00629 void identityMatrix(MultiArrayView<2, T, C> &r)
00630 {
00631     const unsigned int rows = rowCount(r);
00632     vigra_precondition(rows == columnCount(r),
00633        "identityMatrix(): Matrix must be square.");
00634     for(unsigned int i = 0; i < rows; ++i) {
00635         for(unsigned int j = 0; j < rows; ++j)
00636             r(j, i) = NumericTraits<T>::zero();
00637         r(i, i) = NumericTraits<T>::one();
00638     }
00639 }
00640 
00641     /** create n identity matrix of the given size.
00642         Usage:
00643         
00644         \code
00645         vigra::Matrix<double> m = vigra::identityMatrix<double>(size);
00646         \endcode
00647     
00648     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00649     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00650         Namespaces: vigra and vigra::linalg
00651      */ 
00652 template <class T>
00653 TemporaryMatrix<T> identityMatrix(unsigned int size)
00654 {
00655     TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero());
00656     for(unsigned int i = 0; i < size; ++i)
00657         ret(i, i) = NumericTraits<T>::one();
00658     return ret;
00659 }
00660 
00661 template <class T, class C1, class C2>
00662 void diagonalMatrixImpl(MultiArrayView<1, T, C1> const & v, MultiArrayView<2, T, C2> &r)
00663 {
00664     const unsigned int size = v.elementCount();
00665     vigra_precondition(rowCount(r) == size && columnCount(r) == size,
00666         "diagonalMatrix(): result must be a square matrix.");
00667     for(unsigned int i = 0; i < size; ++i)
00668         r(i, i) = v(i);
00669 }
00670 
00671     /** make a diagonal matrix from a vector.
00672         The vector is given as matrix \a v, which must either have a single
00673         row or column. The result is witten into the square matrix \a r.
00674     
00675     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00676     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00677         Namespaces: vigra and vigra::linalg
00678      */ 
00679 template <class T, class C1, class C2>
00680 void diagonalMatrix(MultiArrayView<2, T, C1> const & v, MultiArrayView<2, T, C2> &r)
00681 {
00682     vigra_precondition(rowCount(v) == 1 || columnCount(v) == 1,
00683         "diagonalMatrix(): input must be a vector.");
00684     r.init(NumericTraits<T>::zero());
00685     if(rowCount(v) == 1)
00686         diagonalMatrixImpl(v.bindInner(0), r);
00687     else
00688         diagonalMatrixImpl(v.bindOuter(0), r);
00689 }
00690 
00691     /** create a diagonal matrix from a vector.
00692         The vector is given as matrix \a v, which must either have a single
00693         row or column. The result is returned as a temporary matrix.
00694         Usage:
00695         
00696         \code
00697         vigra::Matrix<double> v(1, len);
00698         v = ...;
00699         
00700         vigra::Matrix<double> m = diagonalMatrix(v);
00701         \endcode
00702     
00703     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00704     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00705         Namespaces: vigra and vigra::linalg
00706      */ 
00707 template <class T, class C>
00708 TemporaryMatrix<T> diagonalMatrix(MultiArrayView<2, T, C> const & v)
00709 {
00710     vigra_precondition(rowCount(v) == 1 || columnCount(v) == 1,
00711         "diagonalMatrix(): input must be a vector.");
00712     unsigned int size = v.elementCount();
00713     TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero());
00714     if(rowCount(v) == 1)
00715         diagonalMatrixImpl(v.bindInner(0), ret);
00716     else
00717         diagonalMatrixImpl(v.bindOuter(0), ret);
00718     return ret;
00719 }
00720 
00721     /** transpose matrix \a v.
00722         The result is written into \a r which must have the correct (i.e.
00723         transposed) shape.
00724     
00725     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00726     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00727         Namespaces: vigra and vigra::linalg
00728      */ 
00729 template <class T, class C1, class C2>
00730 void transpose(const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &r)
00731 {
00732     const unsigned int rows = rowCount(r);
00733     const unsigned int cols = columnCount(r);
00734     vigra_precondition(rows == columnCount(v) && cols == rowCount(v),
00735        "transpose(): arrays must have transposed shapes.");
00736     for(unsigned int i = 0; i < cols; ++i)
00737         for(unsigned int j = 0; j < rows; ++j)
00738             r(j, i) = v(i, j);
00739 }
00740 
00741     /** create the transpose of a matrix \a v.
00742         The result is returned as a temporary matrix.
00743         Usage:
00744         
00745         \code
00746         vigra::Matrix<double> v(rows, cols);
00747         v = ...;
00748         
00749         vigra::Matrix<double> m = transpose(v);
00750         \endcode
00751     
00752     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00753     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00754         Namespaces: vigra and vigra::linalg
00755      */ 
00756 template <class T, class C>
00757 TemporaryMatrix<T> transpose(MultiArrayView<2, T, C> const & v)
00758 {
00759     TemporaryMatrix<T> ret(columnCount(v), rowCount(v));
00760     transpose(v, ret);
00761     return ret;
00762 }
00763 
00764 template <class T>
00765 TemporaryMatrix<T> transpose(TemporaryMatrix<T> const & v)
00766 {
00767     const unsigned int rows = v.rowCount();
00768     const unsigned int cols = v.columnCount();
00769     if(rows == cols)
00770     {
00771         return const_cast<TemporaryMatrix<T> &>(v).transpose();
00772     }
00773     else
00774     {
00775         TemporaryMatrix<T> ret(cols, rows);
00776         transpose(v, ret);
00777         return ret;
00778     }
00779 }
00780 
00781     /** add matrices \a a and \a b.
00782         The result is written into \a r. All three matrices must have the same shape.
00783     
00784     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00785     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00786         Namespace: vigra::linalg
00787      */ 
00788 template <class T, class C1, class C2, class C3>
00789 void add(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
00790               MultiArrayView<2, T, C3> &r)
00791 {
00792     const unsigned int rrows = rowCount(r);
00793     const unsigned int rcols = columnCount(r);
00794     vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) && 
00795                        rrows == rowCount(b) && rcols == columnCount(b),
00796                        "add(): Matrix shapes must agree.");
00797 
00798     for(unsigned int i = 0; i < rcols; ++i) {
00799         for(unsigned int j = 0; j < rrows; ++j) {
00800             r(j, i) = a(j, i) + b(j, i);
00801         }
00802     }
00803 }
00804  
00805     /** add matrices \a a and \a b.
00806         The two matrices must have the same shape.
00807         The result is returned as a temporary matrix. 
00808     
00809     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00810     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00811         Namespace: vigra::linalg
00812      */ 
00813 template <class T, class C1, class C2>
00814 inline TemporaryMatrix<T>
00815 operator+(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
00816 {
00817     return TemporaryMatrix<T>(a) += b;
00818 }
00819 
00820 template <class T, class C>
00821 inline TemporaryMatrix<T>
00822 operator+(const TemporaryMatrix<T> &a, const MultiArrayView<2, T, C> &b)
00823 {
00824     return const_cast<TemporaryMatrix<T> &>(a) += b;
00825 }
00826 
00827 template <class T, class C>
00828 inline TemporaryMatrix<T>
00829 operator+(const MultiArrayView<2, T, C> &a, const TemporaryMatrix<T> &b)
00830 {
00831     return const_cast<TemporaryMatrix<T> &>(b) += a;
00832 }
00833 
00834 template <class T>
00835 inline TemporaryMatrix<T>
00836 operator+(const TemporaryMatrix<T> &a, const TemporaryMatrix<T> &b)
00837 {
00838     return const_cast<TemporaryMatrix<T> &>(a) += b;
00839 }
00840 
00841     /** subtract matrix \a b from \a a.
00842         The result is written into \a r. All three matrices must have the same shape.
00843     
00844     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00845     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00846         Namespace: vigra::linalg
00847      */ 
00848 template <class T, class C1, class C2, class C3>
00849 void sub(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
00850               MultiArrayView<2, T, C3> &r)
00851 {
00852     const unsigned int rrows = rowCount(r);
00853     const unsigned int rcols = columnCount(r);
00854     vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) && 
00855                        rrows == rowCount(b) && rcols == columnCount(b),
00856                        "subtract(): Matrix shapes must agree.");
00857 
00858     for(unsigned int i = 0; i < rcols; ++i) {
00859         for(unsigned int j = 0; j < rrows; ++j) {
00860             r(j, i) = a(j, i) - b(j, i);
00861         }
00862     }
00863 }
00864   
00865     /** subtract matrix \a b from \a a.
00866         The two matrices must have the same shape.
00867         The result is returned as a temporary matrix. 
00868     
00869     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00870     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00871         Namespace: vigra::linalg
00872      */ 
00873 template <class T, class C1, class C2>
00874 inline TemporaryMatrix<T>
00875 operator-(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
00876 {
00877     return TemporaryMatrix<T>(a) -= b;
00878 }
00879 
00880 template <class T, class C>
00881 inline TemporaryMatrix<T>
00882 operator-(const TemporaryMatrix<T> &a, const MultiArrayView<2, T, C> &b)
00883 {
00884     return const_cast<TemporaryMatrix<T> &>(a) -= b;
00885 }
00886 
00887 template <class T, class C>
00888 TemporaryMatrix<T>
00889 operator-(const MultiArrayView<2, T, C> &a, const TemporaryMatrix<T> &b)
00890 {
00891     const unsigned int rows = rowCount(a);
00892     const unsigned int cols = columnCount(a);
00893     vigra_precondition(rows == b.rowCount() && cols == b.columnCount(),
00894        "Matrix::operator-(): Shape mismatch.");
00895     
00896     for(unsigned int i = 0; i < cols; ++i)
00897         for(unsigned int j = 0; j < rows; ++j)
00898             const_cast<TemporaryMatrix<T> &>(b)(j, i) = a(j, i) - b(j, i);
00899     return b;
00900 }
00901 
00902 template <class T>
00903 inline TemporaryMatrix<T>
00904 operator-(const TemporaryMatrix<T> &a, const TemporaryMatrix<T> &b)
00905 {
00906     return const_cast<TemporaryMatrix<T> &>(a) -= b;
00907 }
00908 
00909     /** negate matrix \a a.
00910         The result is returned as a temporary matrix. 
00911     
00912     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00913     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00914         Namespace: vigra::linalg
00915      */ 
00916 template <class T, class C>
00917 inline TemporaryMatrix<T>
00918 operator-(const MultiArrayView<2, T, C> &a)
00919 {
00920     return TemporaryMatrix<T>(a) *= -NumericTraits<T>::one();
00921 }
00922 
00923 template <class T>
00924 inline TemporaryMatrix<T>
00925 operator-(const TemporaryMatrix<T> &a)
00926 {
00927     return const_cast<TemporaryMatrix<T> &>(a) *= -NumericTraits<T>::one();
00928 }
00929 
00930     /** calculate the inner product of two matrices representing vectors. 
00931         That is, matrix \a x must have a single row, and matrix \a y must 
00932         have a single column, and the other dimensions must match.
00933     
00934     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00935     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00936         Namespaces: vigra and vigra::linalg
00937      */ 
00938 template <class T, class C1, class C2>
00939 T dot(const MultiArrayView<2, T, C1> &x, const MultiArrayView<2, T, C2> &y)
00940 {
00941     const unsigned int n = columnCount(x);
00942     vigra_precondition(n == rowCount(y) && 1 == rowCount(x) && 1 == columnCount(y),
00943        "dot(): shape mismatch.");
00944     T ret = NumericTraits<T>::zero();
00945     for(unsigned int i = 0; i < n; ++i)
00946         ret += x(0, i) * y(i, 0);
00947     return ret;
00948 }
00949 
00950     /** calculate the inner product of two vectors. The vector
00951         lenths must match.
00952     
00953     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00954     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00955         Namespaces: vigra and vigra::linalg
00956      */ 
00957 template <class T, class C1, class C2>
00958 T dot(const MultiArrayView<1, T, C1> &x, const MultiArrayView<1, T, C2> &y)
00959 {
00960     const unsigned int n = x.elementCount();
00961     vigra_precondition(n == y.elementCount(),
00962        "dot(): shape mismatch.");
00963     T ret = NumericTraits<T>::zero();
00964     for(unsigned int i = 0; i < n; ++i)
00965         ret += x(i) * y(i);
00966     return ret;
00967 }
00968 
00969     /** calculate the outer product of two matrices representing vectors. 
00970         That is, matrix \a x must have a single column, and matrix \a y must 
00971         have a single row, and the other dimensions must match. The result
00972         is written into \a r.
00973     
00974     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00975     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00976         Namespaces: vigra and vigra::linalg
00977      */ 
00978 template <class T, class C1, class C2, class C3>
00979 void outer(const MultiArrayView<2, T, C1> &x, const MultiArrayView<2, T, C2> &y,
00980       MultiArrayView<2, T, C3> &r)
00981 {
00982     const unsigned int rows = rowCount(r);
00983     const unsigned int cols = columnCount(r);
00984     vigra_precondition(rows == rowCount(x) && cols == columnCount(y) && 
00985                        1 == columnCount(x) && 1 == rowCount(y),
00986        "outer(): shape mismatch.");
00987     for(unsigned int i = 0; i < cols; ++i)
00988         for(unsigned int j = 0; j < rows; ++j)
00989             r(j, i) = x(j, 0) * y(0, i);
00990 }
00991 
00992     /** calculate the outer product of two matrices representing vectors. 
00993         That is, matrix \a x must have a single column, and matrix \a y must 
00994         have a single row, and the other dimensions must match. The result
00995         is returned as a temporary matrix.
00996     
00997     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
00998     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
00999         Namespaces: vigra and vigra::linalg
01000      */ 
01001 template <class T, class C1, class C2>
01002 TemporaryMatrix<T> 
01003 outer(const MultiArrayView<2, T, C1> &x, const MultiArrayView<2, T, C2> &y)
01004 {
01005     const unsigned int rows = rowCount(x);
01006     const unsigned int cols = columnCount(y);
01007     vigra_precondition(1 == columnCount(x) && 1 == rowCount(y),
01008        "outer(): shape mismatch.");
01009     TemporaryMatrix<T> ret(rows, cols);
01010     outer(x, y, ret);
01011     return ret;
01012 }
01013 
01014     /** multiply matrix \a a with scalar \a b.
01015         The result is written into \a r. \a a and \a r must have the same shape.
01016     
01017     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01018     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01019         Namespace: vigra::linalg
01020      */ 
01021 template <class T, class C1, class C2>
01022 void smul(const MultiArrayView<2, T, C1> &a, T b, MultiArrayView<2, T, C2> &r)
01023 {
01024     const unsigned int rows = rowCount(a);
01025     const unsigned int cols = columnCount(a);
01026     vigra_precondition(rows == rowCount(r) && cols == columnCount(r),
01027                        "smul(): Matrix sizes must agree.");
01028     
01029     for(unsigned int i = 0; i < cols; ++i)
01030         for(unsigned int j = 0; j < rows; ++j)
01031             r(j, i) = a(j, i) * b;
01032 }
01033 
01034     /** multiply scalar \a a with matrix \a b.
01035         The result is written into \a r. \a b and \a r must have the same shape.
01036     
01037     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01038     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01039         Namespace: vigra::linalg
01040      */ 
01041 template <class T, class C2, class C3>
01042 void smul(T a, const MultiArrayView<2, T, C2> &b, MultiArrayView<2, T, C3> &r)
01043 {
01044     smul(b, a, r);
01045 }
01046 
01047     /** perform matrix multiplication of matrices \a a and \a b.
01048         The result is written into \a r. The three matrices must have matching shapes.
01049     
01050     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01051     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01052         Namespace: vigra::linalg
01053      */ 
01054 template <class T, class C1, class C2, class C3>
01055 void mmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
01056          MultiArrayView<2, T, C3> &r)
01057 {
01058     const unsigned int rrows = rowCount(r);
01059     const unsigned int rcols = columnCount(r);
01060     const unsigned int acols = columnCount(a);
01061     vigra_precondition(rrows == rowCount(a) && rcols == columnCount(b) && acols == rowCount(b),
01062                        "mmul(): Matrix shapes must agree.");
01063 
01064     for(unsigned int i = 0; i < rcols; ++i) {
01065         for(unsigned int j = 0; j < rrows; ++j) {
01066             r(j, i) = 0.0;
01067             for(unsigned int k = 0; k < acols; ++k) {
01068                 r(j, i) += a(j, k) * b(k, i);
01069             }
01070         }
01071     }
01072 }
01073 
01074     /** perform matrix multiplication of matrices \a a and \a b.
01075         \a a and \a b must have matching shapes.
01076         The result is returned as a temporary matrix. 
01077     
01078     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01079     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01080         Namespace: vigra::linalg
01081      */ 
01082 template <class T, class C1, class C2>
01083 inline TemporaryMatrix<T>
01084 mmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
01085 {
01086     TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
01087     mmul(a, b, ret);
01088     return ret;
01089 }
01090 
01091     /** multiply two matrices \a a and \a b pointwise.
01092         The result is written into \a r. All three matrices must have the same shape.
01093     
01094     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01095     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01096         Namespace: vigra::linalg
01097      */ 
01098 template <class T, class C1, class C2, class C3>
01099 void pmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
01100               MultiArrayView<2, T, C3> &r)
01101 {
01102     const unsigned int rrows = rowCount(r);
01103     const unsigned int rcols = columnCount(r);
01104     vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) && 
01105                        rrows == rowCount(b) && rcols == columnCount(b),
01106                        "pmul(): Matrix shapes must agree.");
01107 
01108     for(unsigned int i = 0; i < rcols; ++i) {
01109         for(unsigned int j = 0; j < rrows; ++j) {
01110             r(j, i) = a(j, i) * b(j, i);
01111         }
01112     }
01113 }
01114 
01115     /** multiply matrices \a a and \a b pointwise.
01116         \a a and \a b must have matching shapes.
01117         The result is returned as a temporary matrix. 
01118     
01119     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01120     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01121         Namespace: vigra::linalg
01122      */ 
01123 template <class T, class C1, class C2>
01124 inline TemporaryMatrix<T>
01125 pmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
01126 {
01127     TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
01128     pmul(a, b, ret);
01129     return ret;
01130 }
01131 
01132     /** multiply matrix \a a with scalar \a b.
01133         The result is returned as a temporary matrix. 
01134     
01135     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01136     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01137         Namespace: vigra::linalg
01138      */ 
01139 template <class T, class C>
01140 inline TemporaryMatrix<T>
01141 operator*(const MultiArrayView<2, T, C> &a, T b)
01142 {
01143     return TemporaryMatrix<T>(a) *= b;
01144 }
01145 
01146 template <class T>
01147 inline TemporaryMatrix<T>
01148 operator*(const TemporaryMatrix<T> &a, T b)
01149 {
01150     return const_cast<TemporaryMatrix<T> &>(a) *= b;
01151 }
01152 
01153     /** multiply scalar \a a with matrix \a b.
01154         The result is returned as a temporary matrix. 
01155     
01156     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01157     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01158         Namespace: vigra::linalg
01159      */ 
01160 template <class T, class C>
01161 inline TemporaryMatrix<T>
01162 operator*(T a, const MultiArrayView<2, T, C> &b)
01163 {
01164     return TemporaryMatrix<T>(b) *= a;
01165 }
01166 
01167 template <class T>
01168 inline TemporaryMatrix<T>
01169 operator*(T a, const TemporaryMatrix<T> &b)
01170 {
01171     return const_cast<TemporaryMatrix<T> &>(b) *= b;
01172 }
01173 
01174     /** perform matrix multiplication of matrices \a a and \a b.
01175         \a a and \a b must have matching shapes.
01176         The result is returned as a temporary matrix. 
01177     
01178     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01179     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01180         Namespace: vigra::linalg
01181      */ 
01182 template <class T, class C1, class C2>
01183 inline TemporaryMatrix<T>
01184 operator*(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
01185 {
01186     TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
01187     mmul(a, b, ret);
01188     return ret;
01189 }
01190 
01191     /** divide matrix \a a by scalar \a b.
01192         The result is written into \a r. \a a and \a r must have the same shape.
01193     
01194     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01195     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01196         Namespace: vigra::linalg
01197      */ 
01198 template <class T, class C1, class C2>
01199 void sdiv(const MultiArrayView<2, T, C1> &a, T b, MultiArrayView<2, T, C2> &r)
01200 {
01201     const unsigned int rows = rowCount(a);
01202     const unsigned int cols = columnCount(a);
01203     vigra_precondition(rows == rowCount(r) && cols == columnCount(r),
01204                        "sdiv(): Matrix sizes must agree.");
01205     
01206     for(unsigned int i = 0; i < cols; ++i)
01207         for(unsigned int j = 0; j < rows; ++j)
01208             r(j, i) = a(j, i) / b;
01209 }
01210 
01211     /** divide two matrices \a a and \a b pointwise.
01212         The result is written into \a r. All three matrices must have the same shape.
01213     
01214     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01215     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01216         Namespace: vigra::linalg
01217      */ 
01218 template <class T, class C1, class C2, class C3>
01219 void pdiv(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
01220               MultiArrayView<2, T, C3> &r)
01221 {
01222     const unsigned int rrows = rowCount(r);
01223     const unsigned int rcols = columnCount(r);
01224     vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) && 
01225                        rrows == rowCount(b) && rcols == columnCount(b),
01226                        "pdiv(): Matrix shapes must agree.");
01227 
01228     for(unsigned int i = 0; i < rcols; ++i) {
01229         for(unsigned int j = 0; j < rrows; ++j) {
01230             r(j, i) = a(j, i) * b(j, i);
01231         }
01232     }
01233 }
01234 
01235     /** divide matrices \a a and \a b pointwise.
01236         \a a and \a b must have matching shapes.
01237         The result is returned as a temporary matrix. 
01238     
01239     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01240     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01241         Namespace: vigra::linalg
01242      */ 
01243 template <class T, class C1, class C2>
01244 inline TemporaryMatrix<T>
01245 pdiv(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
01246 {
01247     TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
01248     pdiv(a, b, ret);
01249     return ret;
01250 }
01251 
01252     /** divide matrix \a a by scalar \a b.
01253         The result is returned as a temporary matrix. 
01254     
01255     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01256     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01257         Namespace: vigra::linalg
01258      */ 
01259 template <class T, class C>
01260 inline TemporaryMatrix<T>
01261 operator/(const MultiArrayView<2, T, C> &a, T b)
01262 {
01263     return TemporaryMatrix<T>(a) /= b;
01264 }
01265 
01266 template <class T>
01267 inline TemporaryMatrix<T>
01268 operator/(const TemporaryMatrix<T> &a, T b)
01269 {
01270     return const_cast<TemporaryMatrix<T> &>(a) /= b;
01271 }
01272 
01273 //@}
01274 
01275 } // namespace linalg
01276 
01277 using linalg::Matrix;
01278 using linalg::identityMatrix;
01279 using linalg::diagonalMatrix;
01280 using linalg::squaredNorm;
01281 using linalg::norm;
01282 using linalg::transpose;
01283 using linalg::dot;
01284 using linalg::outer;
01285 using linalg::rowCount;
01286 using linalg::columnCount;
01287 using linalg::rowVector;
01288 using linalg::columnVector;
01289 using linalg::isSymmetric;
01290 
01291 } // namespace vigra
01292 
01293 namespace std {
01294 
01295 /** \addtogroup LinearAlgebraFunctions Matrix functions
01296  */
01297 //@{
01298  
01299     /** print a matrix \a m to the stream \a s. 
01300     
01301     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01302     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01303         Namespace: std
01304      */ 
01305 template <class T, class C>
01306 ostream &
01307 operator<<(ostream & s, const vigra::MultiArrayView<2, T, C> &m)
01308 {
01309     const unsigned int rows = vigra::linalg::rowCount(m);
01310     const unsigned int cols = vigra::linalg::columnCount(m);
01311     ios::fmtflags flags = 
01312         s.setf(ios::right | ios::fixed, ios::adjustfield | ios::floatfield);
01313     for(unsigned int j = 0; j < rows; ++j) 
01314     {
01315         for(unsigned int i = 0; i < cols; ++i)
01316         {
01317             s << setw(7) << setprecision(4) << m(j, i) << " ";
01318         }
01319         s << endl;
01320     }
01321     s.setf(flags);
01322     return s;
01323 }
01324 
01325 //@}
01326 
01327 } // namespace std
01328 
01329 
01330 
01331 #endif // VIGRA_MATRIX_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.2 (27 Jan 2005)