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