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