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