core/vnl/vnl_matrix.h
Go to the documentation of this file.
00001 // This is core/vnl/vnl_matrix.h
00002 #ifndef vnl_matrix_h_
00003 #define vnl_matrix_h_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief An ordinary mathematical matrix
00010 // \verbatim
00011 //  Modifications
00012 //   Apr 21, 1989 - MBN - Initial design and implementation
00013 //   Jun 22, 1989 - MBN - Removed non-destructive methods
00014 //   Aug 09, 1989 - LGO - Inherit from Generic
00015 //   Aug 20, 1989 - MBN - Changed template usage to reflect new syntax
00016 //   Sep 11, 1989 - MBN - Added conditional exception handling and base class
00017 //   Oct 05, 1989 - LGO - Don't re-allocate data in operator= when same size
00018 //   Oct 19, 1989 - LGO - Add extra parameter to varargs constructor
00019 //   Oct 19, 1989 - MBN - Added optional argument to set_compare method
00020 //   Dec 08, 1989 - LGO - Allocate column data in one chunk
00021 //   Dec 08, 1989 - LGO - Clean-up get and put, add const everywhere.
00022 //   Dec 19, 1989 - LGO - Remove the map and reduce methods
00023 //   Feb 22, 1990 - MBN - Changed size arguments from int to unsigned int
00024 //   Jun 30, 1990 - MJF - Added base class name to constructor initializer
00025 //   Feb 21, 1992 - VDN - New lite version
00026 //   May 05, 1992 - VDN - Use envelope to avoid unnecessary copying
00027 //   Sep 30, 1992 - VDN - Matrix inversion with singular value decomposition
00028 //   Aug 21, 1996 - AWF - set_identity, normalize_rows, scale_row.
00029 //   Sep 30, 1996 - AWF - set_row/column methods. Const-correct data_block().
00030 //   14 Feb 1997  - AWF - get_n_rows, get_n_columns.
00031 //   20 Mar 1997  - PVR - get_row, get_column.
00032 //   24-Oct-2010 - Peter Vanroose - mutators and filling methods now return *this
00033 //   18-Jan-2011 - Peter Vanroose - added methods set_diagonal() & get_diagonal()
00034 // \endverbatim
00035 
00036 #include <vcl_iosfwd.h>
00037 #include <vnl/vnl_tag.h>
00038 #include <vnl/vnl_c_vector.h>
00039 #include <vnl/vnl_config.h>
00040 #ifndef NDEBUG
00041 # if VNL_CONFIG_CHECK_BOUNDS
00042 #  include <vnl/vnl_error.h>
00043 #  include <vcl_cassert.h>
00044 # endif
00045 #else
00046 # undef VNL_CONFIG_CHECK_BOUNDS
00047 # define VNL_CONFIG_CHECK_BOUNDS 0
00048 # undef ERROR_CHECKING
00049 #endif
00050 
00051 export template <class T> class vnl_vector;
00052 export template <class T> class vnl_matrix;
00053 
00054 //--------------------------------------------------------------------------------
00055 
00056 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00057 #define v vnl_vector<T>
00058 #define m vnl_matrix<T>
00059 #endif // DOXYGEN_SHOULD_SKIP_THIS
00060 template <class T> m operator+(T const&, m const&);
00061 template <class T> m operator-(T const&, m const&);
00062 template <class T> m operator*(T const&, m const&);
00063 template <class T> m element_product(m const&, m const&);
00064 template <class T> m element_quotient(m const&, m const&);
00065 template <class T> T dot_product(m const&, m const&);
00066 template <class T> T inner_product(m const&, m const&);
00067 template <class T> T cos_angle(m const&, m const& );
00068 template <class T> vcl_ostream& operator<<(vcl_ostream&, m const&);
00069 template <class T> vcl_istream& operator>>(vcl_istream&, m&);
00070 #undef v
00071 #undef m
00072 
00073 //--------------------------------------------------------------------------------
00074 
00075 enum vnl_matrix_type
00076 {
00077   vnl_matrix_null,
00078   vnl_matrix_identity
00079 };
00080 
00081 //:  An ordinary mathematical matrix
00082 // The vnl_matrix<T> class implements two-dimensional arithmetic
00083 // matrices  for  a user-specified numeric data type. Using the
00084 // parameterized types facility of C++,  it  is  possible,  for
00085 // example, for the user to create a matrix of rational numbers
00086 // by parameterizing the vnl_matrix class over the Rational  class.
00087 // The  only  requirement  for the type is that it supports the
00088 // basic arithmetic operators.
00089 //
00090 // Note: Unlike   the   other   sequence   classes,   the
00091 // vnl_matrix<T>  class is fixed-size. It will not grow once the
00092 // size has been specified to the constructor or changed by the
00093 // assignment  or  multiplication  operators.  The vnl_matrix<T>
00094 // class is row-based with addresses of rows being cached,  and
00095 // elements accessed as m[row][col].
00096 //
00097 // Note: The matrix can, however, be resized using the set_size(nr,nc) function.
00098 //
00099 // Note: Indexing of the matrix is zero-based, so the top-left element is M(0,0).
00100 //
00101 // Note: Inversion of matrix M, and other operations such as solving systems of linear
00102 // equations are handled by the matrix decomposition classes in vnl/algo, such
00103 // as matrix_inverse, svd, qr etc.
00104 //
00105 // Note: Use a vnl_vector<T> with these matrices.
00106 
00107 template<class T>
00108 class vnl_matrix
00109 {
00110  public:
00111   //: Default constructor creates an empty matrix of size 0,0.
00112   vnl_matrix() :
00113     num_rows(0),
00114     num_cols(0),
00115     data(0)
00116   {
00117   }
00118 
00119   //: Construct a matrix of size r rows by c columns
00120   // Contents are unspecified.
00121   // Complexity $O(1)$
00122   vnl_matrix(unsigned r, unsigned c);                           // r rows, c cols.
00123 
00124   //: Construct a matrix of size r rows by c columns, and all elements equal to v0
00125   // Complexity $O(r.c)$
00126   vnl_matrix(unsigned r, unsigned c, T const& v0);              // r rows, c cols, value v0.
00127 
00128   //: Construct a matrix of size r rows by c columns, with a special type
00129   // Contents are specified by t
00130   // Complexity $O(r.c)$
00131   vnl_matrix(unsigned r, unsigned c, vnl_matrix_type t);        // r rows, c cols, special type
00132 
00133   //: Construct a matrix of size r rows by c columns, initialised by an automatic array
00134   // The first n elements, are initialised row-wise, to values.
00135   // Complexity $O(n)$
00136   vnl_matrix(unsigned r, unsigned c, unsigned n, T const values[]);  // use automatic arrays.
00137 
00138   //: Construct a matrix of size r rows by c columns, initialised by a memory block
00139   // The values are initialise row wise from the data.
00140   // Complexity $O(r.c)$
00141   vnl_matrix(T const* data_block, unsigned r, unsigned c);      // fill row-wise.
00142 
00143   //: Copy construct a matrix
00144   // Complexity $O(r.c)$
00145   vnl_matrix(vnl_matrix<T> const&);                             // from another matrix.
00146 
00147 #ifndef VXL_DOXYGEN_SHOULD_SKIP_THIS
00148 // <internal>
00149   // These constructors are here so that operator* etc can take
00150   // advantage of the C++ return value optimization.
00151   vnl_matrix(vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_add); // M + M
00152   vnl_matrix(vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_sub); // M - M
00153   vnl_matrix(vnl_matrix<T> const &, T,                     vnl_tag_mul); // M * s
00154   vnl_matrix(vnl_matrix<T> const &, T,                     vnl_tag_div); // M / s
00155   vnl_matrix(vnl_matrix<T> const &, T,                     vnl_tag_add); // M + s
00156   vnl_matrix(vnl_matrix<T> const &, T,                     vnl_tag_sub); // M - s
00157   vnl_matrix(vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_mul); // M * M
00158   vnl_matrix(vnl_matrix<T> &that, vnl_tag_grab)
00159     : num_rows(that.num_rows), num_cols(that.num_cols), data(that.data)
00160   { that.num_cols=that.num_rows=0; that.data=0; } // "*this" now uses "that"'s data.
00161 // </internal>
00162 #endif
00163 
00164   //: Matrix destructor
00165   ~vnl_matrix();
00166 
00167 // Basic 2D-Array functionality-------------------------------------------
00168 
00169   //: Return number of rows
00170   unsigned rows()    const { return num_rows; }
00171 
00172   //: Return number of columns
00173   // A synonym for cols()
00174   unsigned columns()  const { return num_cols; }
00175 
00176   //: Return number of columns
00177   // A synonym for columns()
00178   unsigned cols()    const { return num_cols; }
00179 
00180   //: Return number of elements
00181   // This equals rows() * cols()
00182   unsigned size()    const { return rows()*cols(); }
00183 
00184   //: set element with boundary checks if error checking is on.
00185   void put(unsigned r, unsigned c, T const&);
00186 
00187   //: get element with boundary checks if error checking is on.
00188   T    get(unsigned r, unsigned c) const;
00189 
00190   //: return pointer to given row
00191   // No boundary checking here.
00192   T       * operator[](unsigned r) { return data[r]; }
00193 
00194   //: return pointer to given row
00195   // No boundary checking here.
00196   T const * operator[](unsigned r) const { return data[r]; }
00197 
00198   //: Access an element for reading or writing
00199   // There are assert style boundary checks - #define NDEBUG to turn them off.
00200   T       & operator()(unsigned r, unsigned c)
00201   {
00202 #if VNL_CONFIG_CHECK_BOUNDS
00203     assert(r<rows());   // Check the row index is valid
00204     assert(c<cols());   // Check the column index is valid
00205 #endif
00206     return this->data[r][c];
00207   }
00208 
00209   //: Access an element for reading
00210   // There are assert style boundary checks - #define NDEBUG to turn them off.
00211   T const & operator()(unsigned r, unsigned c) const
00212   {
00213 #if VNL_CONFIG_CHECK_BOUNDS
00214     assert(r<rows());   // Check the row index is valid
00215     assert(c<cols());   // Check the column index is valid
00216 #endif
00217     return this->data[r][c];
00218   }
00219 
00220 
00221   // ----------------------- Filling and copying -----------------------
00222 
00223   //: Sets all elements of matrix to specified value, and returns "*this".
00224   //  Complexity $O(r.c)$
00225   //  Returning "*this" allows "chaining" two or more operations:
00226   //  e.g., to set a matrix to a column-normalized all-elements-equal matrix, say
00227   //  \code
00228   //     M.fill(1).normalize_columns();
00229   //  \endcode
00230   //  Returning "*this" also allows passing such a matrix as argument
00231   //  to a function f, without having to name the constructed matrix:
00232   //  \code
00233   //     f(vnl_matrix<double>(5,5,1.0).normalize_columns());
00234   //  \endcode
00235   vnl_matrix& fill(T const&);
00236 
00237   //: Sets all diagonal elements of matrix to specified value; returns "*this".
00238   //  Complexity $O(\min(r,c))$
00239   //  Returning "*this" allows "chaining" two or more operations:
00240   //  e.g., to set a 3x3 matrix to [5 0 0][0 10 0][0 0 15], just say
00241   //  \code
00242   //     M.fill_diagonal(5).scale_row(1,2).scale_column(2,3);
00243   //  \endcode
00244   //  Returning "*this" also allows passing a diagonal-filled matrix as argument
00245   //  to a function f, without having to name the constructed matrix:
00246   //  \code
00247   //     f(vnl_matrix<double>(3,3).fill_diagonal(5));
00248   //  \endcode
00249   vnl_matrix& fill_diagonal(T const&);
00250 
00251   //: Sets the diagonal elements of this matrix to the specified list of values.
00252   //  Returning "*this" allows "chaining" two or more operations: see the
00253   //  reasoning (and the examples) in the documentation for method
00254   //  fill_diagonal().
00255   vnl_matrix& set_diagonal(vnl_vector<T> const&);
00256 
00257   //: Fills (laminates) this matrix with the given data, then returns it.
00258   //  We assume that the argument points to a contiguous rows*cols array, stored rowwise.
00259   //  No bounds checking on the array.
00260   //  Returning "*this" allows "chaining" two or more operations:
00261   //  e.g., to fill a square matrix column-wise, fill it rowwise then transpose:
00262   //  \code
00263   //     M.copy_in(array).inplace_transpose();
00264   //  \endcode
00265   //  Returning "*this" also allows passing a filled-in matrix as argument
00266   //  to a function f, without having to name the constructed matrix:
00267   //  \code
00268   //     f(vnl_matrix<double>(3,3).copy_in(array));
00269   //  \endcode
00270   vnl_matrix& copy_in(T const *);
00271 
00272   //: Fills (laminates) this matrix with the given data, then returns it.
00273   // A synonym for copy_in()
00274   vnl_matrix& set(T const *d) { return copy_in(d); }
00275 
00276   //: Fills the given array with this matrix.
00277   //  We assume that the argument points to a contiguous rows*cols array, stored rowwise.
00278   // No bounds checking on the array.
00279   void copy_out(T *) const;
00280 
00281   //: Set all elements to value v
00282   // Complexity $O(r.c)$
00283   vnl_matrix<T>& operator=(T const&v) { fill(v); return *this; }
00284 
00285   //: Copies all elements of rhs matrix into lhs matrix.
00286   // Complexity $O(\min(r,c))$
00287   vnl_matrix<T>& operator=(vnl_matrix<T> const&);
00288 
00289   // ----------------------- Arithmetic --------------------------------
00290   // note that these functions should not pass scalar as a const&.
00291   // Look what would happen to A /= A(0,0).
00292 
00293   //: Add rhs to each element of lhs matrix in situ
00294   vnl_matrix<T>& operator+=(T value);
00295 
00296   //: Subtract rhs from each element of lhs matrix in situ
00297   vnl_matrix<T>& operator-=(T value);
00298 
00299   //: Scalar multiplication in situ of lhs matrix  by rhs
00300   vnl_matrix<T>& operator*=(T value);
00301 
00302   //: Scalar division of lhs matrix  in situ by rhs
00303   vnl_matrix<T>& operator/=(T value);
00304 
00305   //: Add rhs to lhs  matrix in situ
00306   vnl_matrix<T>& operator+=(vnl_matrix<T> const&);
00307   //: Subtract rhs from lhs matrix in situ
00308   vnl_matrix<T>& operator-=(vnl_matrix<T> const&);
00309   //: Multiply lhs matrix in situ by rhs
00310   vnl_matrix<T>& operator*=(vnl_matrix<T> const&rhs) { return *this = (*this) * rhs; }
00311 
00312   //: Negate all elements of matrix
00313   vnl_matrix<T> operator-() const;
00314 
00315 
00316   //: Add rhs to each element of lhs matrix and return result in new matrix
00317   vnl_matrix<T> operator+(T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_add()); }
00318 
00319   //: Subtract rhs from each element of lhs matrix and return result in new matrix
00320   vnl_matrix<T> operator-(T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_sub()); }
00321 
00322   //: Scalar multiplication of lhs matrix by rhs  and return result in new matrix
00323   vnl_matrix<T> operator*(T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_mul()); }
00324 
00325   //: Scalar division of lhs matrix by rhs and return result in new matrix
00326   vnl_matrix<T> operator/(T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_div()); }
00327 
00328   //: Matrix add rhs to lhs matrix and return result in new matrix
00329   vnl_matrix<T> operator+(vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_add()); }
00330   //: Matrix subtract rhs from lhs and return result in new matrix
00331   vnl_matrix<T> operator-(vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_sub()); }
00332   //: Matrix multiply lhs by rhs matrix and return result in new matrix
00333   vnl_matrix<T> operator*(vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_mul()); }
00334 
00335   ////--------------------------- Additions ----------------------------
00336 
00337   //: Make a new matrix by applying function to each element.
00338   vnl_matrix<T> apply(T (*f)(T)) const;
00339 
00340   //: Make a new matrix by applying function to each element.
00341   vnl_matrix<T> apply(T (*f)(T const&)) const;
00342 
00343   //: Return transpose
00344   vnl_matrix<T> transpose() const;
00345 
00346   //: Return conjugate transpose
00347   vnl_matrix<T> conjugate_transpose() const;
00348 
00349   //: Set values of this matrix to those of M, starting at [top,left]
00350   vnl_matrix<T>& update(vnl_matrix<T> const&, unsigned top=0, unsigned left=0);
00351 
00352   //: Set the elements of the i'th column to v[i]  (No bounds checking)
00353   vnl_matrix& set_column(unsigned i, T const * v);
00354 
00355   //: Set the elements of the i'th column to value, then return *this.
00356   vnl_matrix& set_column(unsigned i, T value );
00357 
00358   //: Set j-th column to v, then return *this.
00359   vnl_matrix& set_column(unsigned j, vnl_vector<T> const& v);
00360 
00361   //: Set columns to those in M, starting at starting_column, then return *this.
00362   vnl_matrix& set_columns(unsigned starting_column, vnl_matrix<T> const& M);
00363 
00364   //: Set the elements of the i'th row to v[i]  (No bounds checking)
00365   vnl_matrix& set_row(unsigned i, T const * v);
00366 
00367   //: Set the elements of the i'th row to value, then return *this.
00368   vnl_matrix& set_row(unsigned i, T value );
00369 
00370   //: Set the i-th row
00371   vnl_matrix& set_row(unsigned i, vnl_vector<T> const&);
00372 
00373   //: Extract a sub-matrix of size r x c, starting at (top,left)
00374   //  Thus it contains elements  [top,top+r-1][left,left+c-1]
00375   vnl_matrix<T> extract(unsigned r, unsigned c,
00376                         unsigned top=0, unsigned left=0) const;
00377 
00378   //: Extract a sub-matrix starting at (top,left)
00379   //
00380   //  The output is stored in \a sub_matrix, and it should have the
00381   //  required size on entry.  Thus the result will contain elements
00382   //  [top,top+sub_matrix.rows()-1][left,left+sub_matrix.cols()-1]
00383   void extract ( vnl_matrix<T>& sub_matrix,
00384                  unsigned top=0, unsigned left=0) const;
00385 
00386 
00387   //: Get a vector equal to the given row
00388   vnl_vector<T> get_row(unsigned r) const;
00389 
00390   //: Get a vector equal to the given column
00391   vnl_vector<T> get_column(unsigned c) const;
00392 
00393   //: Get n rows beginning at rowstart
00394   vnl_matrix<T> get_n_rows(unsigned rowstart, unsigned n) const;
00395 
00396   //: Get n columns beginning at colstart
00397   vnl_matrix<T> get_n_columns(unsigned colstart, unsigned n) const;
00398 
00399   //: Return a vector with the content of the (main) diagonal
00400   vnl_vector<T> get_diagonal() const;
00401 
00402   // ==== mutators ====
00403 
00404   //: Sets this matrix to an identity matrix, then returns "*this".
00405   //  Returning "*this" allows e.g. passing an identity matrix as argument to
00406   //  a function f, without having to name the constructed matrix:
00407   //  \code
00408   //     f(vnl_matrix<double>(5,5).set_identity());
00409   //  \endcode
00410   //  Returning "*this" also allows "chaining" two or more operations:
00411   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00412   //  \code
00413   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00414   //  \endcode
00415   //  If the matrix is not square, anyhow set main diagonal to 1, the rest to 0.
00416   vnl_matrix& set_identity();
00417 
00418   //: Transposes this matrix efficiently, and returns it.
00419   //  Returning "*this" allows "chaining" two or more operations:
00420   //  e.g., to fill a square matrix column-wise, fill it rowwise then transpose:
00421   //  \code
00422   //     M.copy_in(array).inplace_transpose();
00423   //  \endcode
00424   vnl_matrix& inplace_transpose();
00425 
00426   //: Reverses the order of rows, and returns "*this".
00427   //  Returning "*this" allows "chaining" two or more operations:
00428   //  e.g., to flip both up-down and left-right, one could just say
00429   //  \code
00430   //     M.flipud().fliplr();
00431   //  \endcode
00432   vnl_matrix& flipud();
00433 
00434   //: Reverses the order of columns, and returns "*this".
00435   //  Returning "*this" allows "chaining" two or more operations:
00436   //  e.g., to flip both up-down and left-right, one could just say
00437   //  \code
00438   //     M.flipud().fliplr();
00439   //  \endcode
00440   vnl_matrix& fliplr();
00441 
00442   //: Normalizes each row so it is a unit vector, and returns "*this".
00443   //  Zero rows are not modified
00444   //  Returning "*this" allows "chaining" two or more operations:
00445   //  e.g., to set a matrix to a row-normalized all-elements-equal matrix, say
00446   //  \code
00447   //     M.fill(1).normalize_rows();
00448   //  \endcode
00449   //  Returning "*this" also allows passing such a matrix as argument
00450   //  to a function f, without having to name the constructed matrix:
00451   //  \code
00452   //     f(vnl_matrix<double>(5,5,1.0).normalize_rows());
00453   //  \endcode
00454   vnl_matrix& normalize_rows();
00455 
00456   //: Normalizes each column so it is a unit vector, and returns "*this".
00457   //  Zero columns are not modified
00458   //  Returning "*this" allows "chaining" two or more operations:
00459   //  e.g., to set a matrix to a column-normalized all-elements-equal matrix, say
00460   //  \code
00461   //     M.fill(1).normalize_columns();
00462   //  \endcode
00463   //  Returning "*this" also allows passing such a matrix as argument
00464   //  to a function f, without having to name the constructed matrix:
00465   //  \code
00466   //     f(vnl_matrix<double>(5,5,1.0).normalize_columns());
00467   //  \endcode
00468   vnl_matrix& normalize_columns();
00469 
00470   //: Scales elements in given row by a factor T, and returns "*this".
00471   //  Returning "*this" allows "chaining" two or more operations:
00472   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00473   //  \code
00474   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00475   //  \endcode
00476   vnl_matrix& scale_row(unsigned row, T value);
00477 
00478   //: Scales elements in given column by a factor T, and returns "*this".
00479   //  Returning "*this" allows "chaining" two or more operations:
00480   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00481   //  \code
00482   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00483   //  \endcode
00484   vnl_matrix& scale_column(unsigned col, T value);
00485 
00486   //: Swap this matrix with that matrix
00487   void swap(vnl_matrix<T> & that);
00488 
00489   //: Type def for norms.
00490   typedef typename vnl_c_vector<T>::abs_t abs_t;
00491 
00492   //: Return sum of absolute values of elements
00493   abs_t array_one_norm() const { return vnl_c_vector<T>::one_norm(begin(), size()); }
00494 
00495   //: Return square root of sum of squared absolute element values
00496   abs_t array_two_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); }
00497 
00498   //: Return largest absolute element value
00499   abs_t array_inf_norm() const { return vnl_c_vector<T>::inf_norm(begin(), size()); }
00500 
00501   //: Return sum of absolute values of elements
00502   abs_t absolute_value_sum() const { return array_one_norm(); }
00503 
00504   //: Return largest absolute value
00505   abs_t absolute_value_max() const { return array_inf_norm(); }
00506 
00507   // $ || M ||_1 := \max_j \sum_i | M_{ij} | $
00508   abs_t operator_one_norm() const;
00509 
00510   // $ || M ||_\inf := \max_i \sum_j | M_{ij} | $
00511   abs_t operator_inf_norm() const;
00512 
00513   //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements)
00514   abs_t frobenius_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); }
00515 
00516   //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements)
00517   abs_t fro_norm() const { return frobenius_norm(); }
00518 
00519   //: Return RMS of all elements
00520   abs_t rms() const { return vnl_c_vector<T>::rms_norm(begin(), size()); }
00521 
00522   //: Return minimum value of elements
00523   T min_value() const { return vnl_c_vector<T>::min_value(begin(), size()); }
00524 
00525   //: Return maximum value of elements
00526   T max_value() const { return vnl_c_vector<T>::max_value(begin(), size()); }
00527 
00528   //: Return location of minimum value of elements
00529   unsigned arg_min() const { return vnl_c_vector<T>::arg_min(begin(), size()); }
00530 
00531   //: Return location of maximum value of elements
00532   unsigned arg_max() const { return vnl_c_vector<T>::arg_max(begin(), size()); }
00533 
00534   //: Return mean of all matrix elements
00535   T mean() const { return vnl_c_vector<T>::mean(begin(), size()); }
00536 
00537   // predicates
00538 
00539   //: Return true iff the size is zero.
00540   bool empty() const { return !data || !num_rows || !num_cols; }
00541 
00542   //:  Return true if all elements equal to identity.
00543   bool is_identity() const;
00544 
00545   //:  Return true if all elements equal to identity, within given tolerance
00546   bool is_identity(double tol) const;
00547 
00548   //: Return true if all elements equal to zero.
00549   bool is_zero() const;
00550 
00551   //: Return true if all elements equal to zero, within given tolerance
00552   bool is_zero(double tol) const;
00553 
00554   //:  Return true if all elements of both matrices are equal, within given tolerance
00555   bool is_equal(vnl_matrix<T> const& rhs, double tol) const;
00556 
00557   //: Return true if finite
00558   bool is_finite() const;
00559 
00560   //: Return true if matrix contains NaNs
00561   bool has_nans() const;
00562 
00563   //: abort if size is not as expected
00564   // This function does or tests nothing if NDEBUG is defined
00565   void assert_size(unsigned r, unsigned c) const
00566   {
00567 #ifndef NDEBUG
00568     assert_size_internal(r, c);
00569 #endif
00570   }
00571   //: abort if matrix contains any INFs or NANs.
00572   // This function does or tests nothing if NDEBUG is defined
00573   void assert_finite() const
00574   {
00575 #ifndef NDEBUG
00576     assert_finite_internal();
00577 #endif
00578   }
00579 
00580   ////----------------------- Input/Output ----------------------------
00581 
00582   //: Read a vnl_matrix from an ascii vcl_istream, automatically determining file size if the input matrix has zero size.
00583   static vnl_matrix<T> read(vcl_istream& s);
00584 
00585   // : Read a vnl_matrix from an ascii vcl_istream, automatically determining file size if the input matrix has zero size.
00586   bool read_ascii(vcl_istream& s);
00587 
00588   //--------------------------------------------------------------------------------
00589 
00590   //: Access the contiguous block storing the elements in the matrix row-wise. O(1).
00591   // 1d array, row-major order.
00592   T const* data_block() const { return data[0]; }
00593 
00594   //: Access the contiguous block storing the elements in the matrix row-wise. O(1).
00595   // 1d array, row-major order.
00596   T      * data_block() { return data[0]; }
00597 
00598   //: Access the 2D array, so that elements can be accessed with array[row][col] directly.
00599   //  2d array, [row][column].
00600   T const* const* data_array() const { return data; }
00601 
00602   //: Access the 2D array, so that elements can be accessed with array[row][col] directly.
00603   //  2d array, [row][column].
00604   T      *      * data_array() { return data; }
00605 
00606   typedef T element_type;
00607 
00608   //: Iterators
00609   typedef T       *iterator;
00610   //: Iterator pointing to start of data
00611   iterator       begin() { return data?data[0]:0; }
00612   //: Iterator pointing to element beyond end of data
00613   iterator       end() { return data?data[0]+num_rows*num_cols:0; }
00614 
00615   //: Const iterators
00616   typedef T const *const_iterator;
00617   //: Iterator pointing to start of data
00618   const_iterator begin() const { return data?data[0]:0; }
00619   //: Iterator pointing to element beyond end of data
00620   const_iterator end() const { return data?data[0]+num_rows*num_cols:0; }
00621 
00622   //: Return a reference to this.
00623   // Useful in code which would prefer not to know if its argument
00624   // is a matrix, matrix_ref or a matrix_fixed.  Note that it doesn't
00625   // return a matrix_ref, so it's only useful in templates or macros.
00626   vnl_matrix<T> const& as_ref() const { return *this; }
00627 
00628   //: Return a reference to this.
00629   vnl_matrix<T>&       as_ref()       { return *this; }
00630 
00631   //--------------------------------------------------------------------------------
00632 
00633   //: Return true if *this == rhs
00634   bool operator_eq(vnl_matrix<T> const & rhs) const;
00635 
00636   //: Equality operator
00637   bool operator==(vnl_matrix<T> const &that) const { return  this->operator_eq(that); }
00638 
00639   //: Inequality operator
00640   bool operator!=(vnl_matrix<T> const &that) const { return !this->operator_eq(that); }
00641 
00642   //: Print matrix to os in some hopefully sensible format
00643   void print(vcl_ostream& os) const;
00644 
00645   //: Make the matrix as if it had been default-constructed.
00646   void clear();
00647 
00648   //: Resize to r rows by c columns. Old data lost.
00649   // Returns true if size changed.
00650   bool set_size(unsigned r, unsigned c);
00651 
00652 //--------------------------------------------------------------------------------
00653 
00654  protected:
00655   unsigned num_rows;   // Number of rows
00656   unsigned num_cols;   // Number of columns
00657   T** data;            // Pointer to the vnl_matrix
00658 
00659 #if VCL_HAS_SLICED_DESTRUCTOR_BUG
00660   // Since this bug exists, we need a flag that can be set during
00661   // construction to tell our destructor whether we own data.
00662   char vnl_matrix_own_data;
00663 #endif
00664 
00665   void assert_size_internal(unsigned r, unsigned c) const;
00666   void assert_finite_internal() const;
00667 
00668   //: Delete data
00669   void destroy();
00670 
00671 #if VCL_NEED_FRIEND_FOR_TEMPLATE_OVERLOAD
00672 # define v vnl_vector<T>
00673 # define m vnl_matrix<T>
00674   friend m operator+         VCL_NULL_TMPL_ARGS (T const&, m const&);
00675   friend m operator-         VCL_NULL_TMPL_ARGS (T const&, m const&);
00676   friend m operator*         VCL_NULL_TMPL_ARGS (T const&, m const&);
00677   friend m element_product   VCL_NULL_TMPL_ARGS (m const&, m const&);
00678   friend m element_quotient  VCL_NULL_TMPL_ARGS (m const&, m const&);
00679   friend T dot_product       VCL_NULL_TMPL_ARGS (m const&, m const&);
00680   friend T inner_product     VCL_NULL_TMPL_ARGS (m const&, m const&);
00681   friend T cos_angle         VCL_NULL_TMPL_ARGS (m const&, m const&);
00682   friend vcl_ostream& operator<< VCL_NULL_TMPL_ARGS (vcl_ostream&, m const&);
00683   friend vcl_istream& operator>> VCL_NULL_TMPL_ARGS (vcl_istream&, m&);
00684 # undef v
00685 # undef m
00686 #endif
00687 
00688   // inline function template instantiation hack for gcc 2.97 -- fsm
00689   static void inline_function_tickler();
00690 };
00691 
00692 
00693 // Definitions of inline functions.
00694 
00695 
00696 //: Returns the value of the element at specified row and column. O(1).
00697 // Checks for valid range of indices.
00698 
00699 template<class T>
00700 inline T vnl_matrix<T>::get(unsigned row, unsigned column) const
00701 {
00702 #ifdef ERROR_CHECKING
00703   if (row >= this->num_rows)                   // If invalid size specified
00704     vnl_error_matrix_row_index("get", row);    // Raise exception
00705   if (column >= this->num_cols)                // If invalid size specified
00706     vnl_error_matrix_col_index("get", column); // Raise exception
00707 #endif
00708   return this->data[row][column];
00709 }
00710 
00711 //: Puts value into element at specified row and column. O(1).
00712 // Checks for valid range of indices.
00713 
00714 template<class T>
00715 inline void vnl_matrix<T>::put(unsigned row, unsigned column, T const& value)
00716 {
00717 #ifdef ERROR_CHECKING
00718   if (row >= this->num_rows)                   // If invalid size specified
00719     vnl_error_matrix_row_index("put", row);    // Raise exception
00720   if (column >= this->num_cols)                // If invalid size specified
00721     vnl_error_matrix_col_index("put", column); // Raise exception
00722 #endif
00723   this->data[row][column] = value;             // Assign data value
00724 }
00725 
00726 
00727 // non-member arithmetical operators.
00728 
00729 //:
00730 // \relatesalso vnl_matrix
00731 template<class T>
00732 inline vnl_matrix<T> operator*(T const& value, vnl_matrix<T> const& m)
00733 {
00734   return vnl_matrix<T>(m, value, vnl_tag_mul());
00735 }
00736 
00737 //:
00738 // \relatesalso vnl_matrix
00739 template<class T>
00740 inline vnl_matrix<T> operator+(T const& value, vnl_matrix<T> const& m)
00741 {
00742   return vnl_matrix<T>(m, value, vnl_tag_add());
00743 }
00744 
00745 //: Swap two matrices
00746 // \relatesalso vnl_matrix
00747 template<class T>
00748 inline void swap(vnl_matrix<T> &A, vnl_matrix<T> &B) { A.swap(B); }
00749 
00750 
00751 #endif // vnl_matrix_h_