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