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