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

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