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