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

| vigra/matrix.hxx |
--- | --- | ---

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)