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

details vigra/matrix.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*        Copyright 2004 by Gunnar Kedenburg and Ullrich Koethe         */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    ( Version 1.5.0, Dec 07 2006 )                                    */
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 "multi_array.hxx"
00046 #include "mathutil.hxx"
00047 #include "numerictraits.hxx"
00048 
00049 
00050 namespace vigra
00051 {
00052 
00053 namespace linalg
00054 {
00055 
00056 template <class T, class C>
00057 inline std::size_t rowCount(const MultiArrayView<2, T, C> &x);
00058 
00059 template <class T, class C>
00060 inline std::size_t 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(std::size_t rows, std::size_t 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(std::size_t rows, std::size_t 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(std::size_t rows, std::size_t 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(std::size_t 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(std::size_t d) const
00297     {
00298         return vigra::linalg::columnVector(*this, d);
00299     }
00300 
00301         /** number of rows (height) of the matrix.
00302         */
00303     std::size_t rowCount() const
00304     {
00305         return this->m_shape[0];
00306     }
00307 
00308         /** number of columns (width) of the matrix.
00309         */
00310     std::size_t columnCount() const
00311     {
00312         return this->m_shape[1];
00313     }
00314 
00315         /** number of elements (width*height) of the matrix.
00316         */
00317     std::size_t 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()(std::size_t row, std::size_t 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()(std::size_t row, std::size_t 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 std::size_t cols = columnCount();
00386     vigra_precondition(cols == rowCount(),
00387         "Matrix::transpose(): in-place transposition requires square matrix.");
00388     for(std::size_t i = 0; i < cols; ++i)
00389         for(std::size_t 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 std::size_t rows = rowCount();
00399     const std::size_t cols = columnCount();
00400     vigra_precondition(rows == vigra::linalg::rowCount(other) && cols == vigra::linalg::columnCount(other),
00401        "Matrix::operator+=(): Shape mismatch.");
00402 
00403     for(std::size_t i = 0; i < cols; ++i)
00404         for(std::size_t 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 std::size_t rows = rowCount();
00414     const std::size_t cols = columnCount();
00415     vigra_precondition(rows == vigra::linalg::rowCount(other) && cols == vigra::linalg::columnCount(other),
00416        "Matrix::operator-=(): Shape mismatch.");
00417 
00418     for(std::size_t i = 0; i < cols; ++i)
00419         for(std::size_t 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 std::size_t rows = rowCount();
00428     const std::size_t cols = columnCount();
00429 
00430     for(std::size_t i = 0; i < cols; ++i)
00431         for(std::size_t 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 std::size_t rows = rowCount();
00440     const std::size_t cols = columnCount();
00441 
00442     for(std::size_t i = 0; i < cols; ++i)
00443         for(std::size_t 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(std::size_t rows, std::size_t columns)
00470     : BaseType(rows, columns, ALLOC())
00471     {}
00472 
00473     TemporaryMatrix(std::size_t rows, std::size_t 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 std::size_t 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 std::size_t 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 std::size_t size = rowCount(m);
00591     if(size != columnCount(m))
00592         return false;
00593 
00594     for(std::size_t i = 0; i < size; ++i)
00595         for(std::size_t 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 Frobenius norm of a matrix.
00614         Equal to the root of 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 std::size_t rows = rowCount(r);
00635     vigra_precondition(rows == columnCount(r),
00636        "identityMatrix(): Matrix must be square.");
00637     for(std::size_t i = 0; i < rows; ++i) {
00638         for(std::size_t 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(std::size_t size)
00657 {
00658     TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero());
00659     for(std::size_t 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 std::size_t size = v.elementCount();
00668     vigra_precondition(rowCount(r) == size && columnCount(r) == size,
00669         "diagonalMatrix(): result must be a square matrix.");
00670     for(std::size_t 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     std::size_t 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 std::size_t rows = rowCount(r);
00736     const std::size_t cols = columnCount(r);
00737     vigra_precondition(rows == columnCount(v) && cols == rowCount(v),
00738        "transpose(): arrays must have transposed shapes.");
00739     for(std::size_t i = 0; i < cols; ++i)
00740         for(std::size_t 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 std::size_t rows = v.rowCount();
00771     const std::size_t 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 std::size_t rrows = rowCount(r);
00796     const std::size_t 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(std::size_t i = 0; i < rcols; ++i) {
00802         for(std::size_t 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 std::size_t rrows = rowCount(r);
00856     const std::size_t 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(std::size_t i = 0; i < rcols; ++i) {
00862         for(std::size_t 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 std::size_t rows = rowCount(a);
00895     const std::size_t cols = columnCount(a);
00896     vigra_precondition(rows == b.rowCount() && cols == b.columnCount(),
00897        "Matrix::operator-(): Shape mismatch.");
00898 
00899     for(std::size_t i = 0; i < cols; ++i)
00900         for(std::size_t 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 std::size_t 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(std::size_t 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         lengths 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 std::size_t n = x.elementCount();
00964     vigra_precondition(n == y.elementCount(),
00965        "dot(): shape mismatch.");
00966     T ret = NumericTraits<T>::zero();
00967     for(std::size_t 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 std::size_t rows = rowCount(r);
01039     const std::size_t 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(std::size_t i = 0; i < cols; ++i)
01044         for(std::size_t 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 std::size_t rows = rowCount(x);
01062     const std::size_t 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     /** calculate the outer product of a matrix (representing a vector) with itself.
01071         The result is returned as a temporary matrix.
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         Namespaces: vigra and vigra::linalg
01076      */
01077 template <class T, class C1>
01078 TemporaryMatrix<T>
01079 outer(const MultiArrayView<2, T, C1> &x)
01080 {
01081     const std::size_t rows = rowCount(x);
01082     const std::size_t cols = columnCount(x);
01083     vigra_precondition(rows == 1 || cols == 1,
01084        "outer(): matrix does not represent a vector.");
01085     const std::size_t size = std::max(rows, cols);
01086     TemporaryMatrix<T> ret(size, size);
01087 
01088     if(rows == 1)
01089     {
01090         for(std::size_t i = 0; i < size; ++i)
01091             for(std::size_t j = 0; j < size; ++j)
01092                 ret(j, i) = x(0, j) * x(0, i);
01093     }
01094     else
01095     {
01096         for(std::size_t i = 0; i < size; ++i)
01097             for(std::size_t j = 0; j < size; ++j)
01098                 ret(j, i) = x(j, 0) * x(i, 0);
01099     }
01100     return ret;
01101 }
01102 
01103     /** multiply matrix \a a with scalar \a b.
01104         The result is written into \a r. \a a and \a r must have the same shape.
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>
01111 void smul(const MultiArrayView<2, T, C1> &a, T b, MultiArrayView<2, T, C2> &r)
01112 {
01113     const std::size_t rows = rowCount(a);
01114     const std::size_t cols = columnCount(a);
01115     vigra_precondition(rows == rowCount(r) && cols == columnCount(r),
01116                        "smul(): Matrix sizes must agree.");
01117 
01118     for(std::size_t i = 0; i < cols; ++i)
01119         for(std::size_t j = 0; j < rows; ++j)
01120             r(j, i) = a(j, i) * b;
01121 }
01122 
01123     /** multiply scalar \a a with matrix \a b.
01124         The result is written into \a r. \a b and \a r must have the same shape.
01125 
01126     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01127     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01128         Namespace: vigra::linalg
01129      */
01130 template <class T, class C2, class C3>
01131 void smul(T a, const MultiArrayView<2, T, C2> &b, MultiArrayView<2, T, C3> &r)
01132 {
01133     smul(b, a, r);
01134 }
01135 
01136     /** perform matrix multiplication of matrices \a a and \a b.
01137         The result is written into \a r. The three matrices must have matching shapes.
01138 
01139     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01140     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01141         Namespace: vigra::linalg
01142      */
01143 template <class T, class C1, class C2, class C3>
01144 void mmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
01145          MultiArrayView<2, T, C3> &r)
01146 {
01147     const std::size_t rrows = rowCount(r);
01148     const std::size_t rcols = columnCount(r);
01149     const std::size_t acols = columnCount(a);
01150     vigra_precondition(rrows == rowCount(a) && rcols == columnCount(b) && acols == rowCount(b),
01151                        "mmul(): Matrix shapes must agree.");
01152 
01153     for(std::size_t i = 0; i < rcols; ++i) {
01154         for(std::size_t j = 0; j < rrows; ++j) {
01155             r(j, i) = 0.0;
01156             for(std::size_t k = 0; k < acols; ++k) {
01157                 r(j, i) += a(j, k) * b(k, i);
01158             }
01159         }
01160     }
01161 }
01162 
01163     /** perform matrix multiplication of matrices \a a and \a b.
01164         \a a and \a b must have matching shapes.
01165         The result is returned as a temporary matrix.
01166 
01167     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01168     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01169         Namespace: vigra::linalg
01170      */
01171 template <class T, class C1, class C2>
01172 inline TemporaryMatrix<T>
01173 mmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
01174 {
01175     TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
01176     mmul(a, b, ret);
01177     return ret;
01178 }
01179 
01180     /** multiply two matrices \a a and \a b pointwise.
01181         The result is written into \a r. All three matrices must have the same shape.
01182 
01183     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01184     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01185         Namespace: vigra::linalg
01186      */
01187 template <class T, class C1, class C2, class C3>
01188 void pmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
01189               MultiArrayView<2, T, C3> &r)
01190 {
01191     const std::size_t rrows = rowCount(r);
01192     const std::size_t rcols = columnCount(r);
01193     vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) &&
01194                        rrows == rowCount(b) && rcols == columnCount(b),
01195                        "pmul(): Matrix shapes must agree.");
01196 
01197     for(std::size_t i = 0; i < rcols; ++i) {
01198         for(std::size_t j = 0; j < rrows; ++j) {
01199             r(j, i) = a(j, i) * b(j, i);
01200         }
01201     }
01202 }
01203 
01204     /** multiply matrices \a a and \a b pointwise.
01205         \a a and \a b must have matching shapes.
01206         The result is returned as a temporary matrix.
01207 
01208     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01209     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01210         Namespace: vigra::linalg
01211      */
01212 template <class T, class C1, class C2>
01213 inline TemporaryMatrix<T>
01214 pmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
01215 {
01216     TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
01217     pmul(a, b, ret);
01218     return ret;
01219 }
01220 
01221     /** multiply matrix \a a with scalar \a b.
01222         The result is returned as a temporary matrix.
01223 
01224     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01225     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01226         Namespace: vigra::linalg
01227      */
01228 template <class T, class C>
01229 inline TemporaryMatrix<T>
01230 operator*(const MultiArrayView<2, T, C> &a, T b)
01231 {
01232     return TemporaryMatrix<T>(a) *= b;
01233 }
01234 
01235 template <class T>
01236 inline TemporaryMatrix<T>
01237 operator*(const TemporaryMatrix<T> &a, T b)
01238 {
01239     return const_cast<TemporaryMatrix<T> &>(a) *= b;
01240 }
01241 
01242     /** multiply scalar \a a with matrix \a b.
01243         The result is returned as a temporary matrix.
01244 
01245     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01246     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01247         Namespace: vigra::linalg
01248      */
01249 template <class T, class C>
01250 inline TemporaryMatrix<T>
01251 operator*(T a, const MultiArrayView<2, T, C> &b)
01252 {
01253     return TemporaryMatrix<T>(b) *= a;
01254 }
01255 
01256 template <class T>
01257 inline TemporaryMatrix<T>
01258 operator*(T a, const TemporaryMatrix<T> &b)
01259 {
01260     return const_cast<TemporaryMatrix<T> &>(b) *= b;
01261 }
01262 
01263     /** multiply matrix \a a with TinyVector \a b.
01264         \a a must be of size <tt>N x N</tt>. Vector \a b and the result
01265         vector are interpreted as column vectors.
01266 
01267     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01268     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01269         Namespace: vigra::linalg
01270      */
01271 template <class T, class A, int N, class DATA, class DERIVED>
01272 TinyVector<T, N>
01273 operator*(const Matrix<T, A> &a, const TinyVectorBase<T, N, DATA, DERIVED> &b)
01274 {
01275     vigra_precondition(N == rowCount(a) && N == columnCount(a),
01276          "operator*(Matrix, TinyVector): Shape mismatch.");
01277 
01278     TinyVector<T, N> res = TinyVectorView<T, N>(&a(0,0)) * b[0];
01279     for(std::size_t i = 1; i < N; ++i)
01280         res += TinyVectorView<T, N>(&a(0,i)) * b[i];
01281     return res;
01282 }
01283 
01284     /** multiply TinyVector \a a with matrix \a b.
01285         \a b must be of size <tt>N x N</tt>. Vector \a a and the result
01286         vector are interpreted as row vectors.
01287 
01288     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01289     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01290         Namespace: vigra::linalg
01291      */
01292 template <class T, int N, class DATA, class DERIVED, class A>
01293 TinyVector<T, N>
01294 operator*(const TinyVectorBase<T, N, DATA, DERIVED> &a, const Matrix<T, A> &b)
01295 {
01296     vigra_precondition(N == rowCount(b) && N == columnCount(b),
01297          "operator*(TinyVector, Matrix): Shape mismatch.");
01298 
01299     TinyVector<T, N> res;
01300     for(std::size_t i = 0; i < N; ++i)
01301         res[i] = dot(a, TinyVectorView<T, N>(&b(0,i)));
01302     return res;
01303 }
01304 
01305     /** perform matrix multiplication of matrices \a a and \a b.
01306         \a a and \a b must have matching shapes.
01307         The result is returned as a temporary matrix.
01308 
01309     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01310     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01311         Namespace: vigra::linalg
01312      */
01313 template <class T, class C1, class C2>
01314 inline TemporaryMatrix<T>
01315 operator*(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
01316 {
01317     TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
01318     mmul(a, b, ret);
01319     return ret;
01320 }
01321 
01322     /** divide matrix \a a by scalar \a b.
01323         The result is written into \a r. \a a and \a r must have the same shape.
01324 
01325     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01326     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01327         Namespace: vigra::linalg
01328      */
01329 template <class T, class C1, class C2>
01330 void sdiv(const MultiArrayView<2, T, C1> &a, T b, MultiArrayView<2, T, C2> &r)
01331 {
01332     const std::size_t rows = rowCount(a);
01333     const std::size_t cols = columnCount(a);
01334     vigra_precondition(rows == rowCount(r) && cols == columnCount(r),
01335                        "sdiv(): Matrix sizes must agree.");
01336 
01337     for(std::size_t i = 0; i < cols; ++i)
01338         for(std::size_t j = 0; j < rows; ++j)
01339             r(j, i) = a(j, i) / b;
01340 }
01341 
01342     /** divide two matrices \a a and \a b pointwise.
01343         The result is written into \a r. All three matrices must have the same shape.
01344 
01345     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01346     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01347         Namespace: vigra::linalg
01348      */
01349 template <class T, class C1, class C2, class C3>
01350 void pdiv(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
01351               MultiArrayView<2, T, C3> &r)
01352 {
01353     const std::size_t rrows = rowCount(r);
01354     const std::size_t rcols = columnCount(r);
01355     vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) &&
01356                        rrows == rowCount(b) && rcols == columnCount(b),
01357                        "pdiv(): Matrix shapes must agree.");
01358 
01359     for(std::size_t i = 0; i < rcols; ++i) {
01360         for(std::size_t j = 0; j < rrows; ++j) {
01361             r(j, i) = a(j, i) * b(j, i);
01362         }
01363     }
01364 }
01365 
01366     /** divide matrices \a a and \a b pointwise.
01367         \a a and \a b must have matching shapes.
01368         The result is returned as a temporary matrix.
01369 
01370     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01371     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01372         Namespace: vigra::linalg
01373      */
01374 template <class T, class C1, class C2>
01375 inline TemporaryMatrix<T>
01376 pdiv(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b)
01377 {
01378     TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
01379     pdiv(a, b, ret);
01380     return ret;
01381 }
01382 
01383     /** divide matrix \a a by scalar \a b.
01384         The result is returned as a temporary matrix.
01385 
01386     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01387     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01388         Namespace: vigra::linalg
01389      */
01390 template <class T, class C>
01391 inline TemporaryMatrix<T>
01392 operator/(const MultiArrayView<2, T, C> &a, T b)
01393 {
01394     return TemporaryMatrix<T>(a) /= b;
01395 }
01396 
01397 template <class T>
01398 inline TemporaryMatrix<T>
01399 operator/(const TemporaryMatrix<T> &a, T b)
01400 {
01401     return const_cast<TemporaryMatrix<T> &>(a) /= b;
01402 }
01403 
01404 //@}
01405 
01406 } // namespace linalg
01407 
01408 using linalg::RowMajor;
01409 using linalg::ColumnMajor;
01410 using linalg::Matrix;
01411 using linalg::identityMatrix;
01412 using linalg::diagonalMatrix;
01413 using linalg::transpose;
01414 using linalg::dot;
01415 using linalg::cross;
01416 using linalg::outer;
01417 using linalg::rowCount;
01418 using linalg::columnCount;
01419 using linalg::rowVector;
01420 using linalg::columnVector;
01421 using linalg::isSymmetric;
01422 
01423 /********************************************************/
01424 /*                                                      */
01425 /*                       NormTraits                     */
01426 /*                                                      */
01427 /********************************************************/
01428 
01429 template <class T, class ALLOC>
01430 struct NormTraits<linalg::Matrix<T, ALLOC> >
01431 {
01432     typedef linalg::Matrix<T, ALLOC> Type;
01433     typedef typename Type::SquaredNormType SquaredNormType;
01434     typedef typename Type::NormType NormType;
01435 };
01436 
01437 template <class T, class ALLOC>
01438 struct NormTraits<linalg::TemporaryMatrix<T, ALLOC> >
01439 {
01440     typedef linalg::TemporaryMatrix<T, ALLOC> Type;
01441     typedef typename Type::SquaredNormType SquaredNormType;
01442     typedef typename Type::NormType NormType;
01443 };
01444 
01445 /** \addtogroup LinearAlgebraFunctions Matrix functions
01446  */
01447 //@{
01448 
01449     /** print a matrix \a m to the stream \a s.
01450 
01451     <b>\#include</b> "<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>" or<br>
01452     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
01453         Namespace: std
01454      */
01455 template <class T, class C>
01456 std::ostream &
01457 operator<<(std::ostream & s, const vigra::MultiArrayView<2, T, C> &m)
01458 {
01459     const std::size_t rows = vigra::linalg::rowCount(m);
01460     const std::size_t cols = vigra::linalg::columnCount(m);
01461     std::ios::fmtflags flags =
01462         s.setf(std::ios::right | std::ios::fixed, std::ios::adjustfield | std::ios::floatfield);
01463     for(std::size_t j = 0; j < rows; ++j)
01464     {
01465         for(std::size_t i = 0; i < cols; ++i)
01466         {
01467             s << std::setw(7) << std::setprecision(4) << m(j, i) << " ";
01468         }
01469         s << std::endl;
01470     }
01471     s.setf(flags);
01472     return s;
01473 }
01474 
01475 //@}
01476 
01477 }  // namespace vigra
01478 
01479 
01480 
01481 #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.5.0 (7 Dec 2006)