core/vnl/vnl_matrix_fixed_ref.txx
Go to the documentation of this file.
00001 // This is core/vnl/vnl_matrix_fixed_ref.txx
00002 #ifndef vnl_matrix_fixed_ref_txx_
00003 #define vnl_matrix_fixed_ref_txx_
00004 
00005 #include "vnl_matrix_fixed_ref.h"
00006 //:
00007 // \file
00008 
00009 #include <vcl_cmath.h>
00010 #include <vcl_iostream.h>
00011 #include <vcl_cstdlib.h> // for abort
00012 #include <vcl_cassert.h>
00013 
00014 #include <vnl/vnl_error.h>
00015 #include <vnl/vnl_math.h>
00016 
00017 #if 0 // commented out
00018 template<class T, unsigned nrows, unsigned ncols>
00019 void
00020 vnl_matrix_fixed_ref_const<T,nrows,ncols>::add( const T* a, const T* b, T* r )
00021 {
00022   unsigned int count = nrows*ncols;
00023   while ( count-- )
00024     *(r++) = *(a++) + *(b++);
00025 }
00026 
00027 
00028 template<class T, unsigned nrows, unsigned ncols>
00029 void
00030 vnl_matrix_fixed_ref_const<T,nrows,ncols>::add( const T* a, T b, T* r )
00031 {
00032   unsigned int count = nrows*ncols;
00033   while ( count-- )
00034     *(r++) = *(a++) + b;
00035 }
00036 
00037 template<class T, unsigned nrows, unsigned ncols>
00038 void
00039 vnl_matrix_fixed_ref_const<T,nrows,ncols>::sub( const T* a, const T* b, T* r )
00040 {
00041   unsigned int count = nrows*ncols;
00042   while ( count-- )
00043     *(r++) = *(a++) - *(b++);
00044 }
00045 
00046 template<class T, unsigned nrows, unsigned ncols>
00047 void
00048 vnl_matrix_fixed_ref_const<T,nrows,ncols>::sub( const T* a, T b, T* r )
00049 {
00050   unsigned int count = nrows*ncols;
00051   while ( count-- )
00052     *(r++) = *(a++) - b;
00053 }
00054 
00055 template<class T, unsigned nrows, unsigned ncols>
00056 void
00057 vnl_matrix_fixed_ref_const<T,nrows,ncols>::sub( T a, const T* b, T* r )
00058 {
00059   unsigned int count = nrows*ncols;
00060   while ( count-- )
00061     *(r++) = a - *(b++);
00062 }
00063 
00064 template<class T, unsigned nrows, unsigned ncols>
00065 void
00066 vnl_matrix_fixed_ref_const<T,nrows,ncols>::mul( const T* a, const T* b, T* r )
00067 {
00068   unsigned int count = nrows*ncols;
00069   while ( count-- )
00070     *(r++) = *(a++) * *(b++);
00071 }
00072 
00073 template<class T, unsigned nrows, unsigned ncols>
00074 void
00075 vnl_matrix_fixed_ref_const<T,nrows,ncols>::mul( const T* a, T b, T* r )
00076 {
00077   unsigned int count = nrows*ncols;
00078   while ( count-- )
00079     *(r++) = *(a++) * b;
00080 }
00081 
00082 template<class T, unsigned nrows, unsigned ncols>
00083 void
00084 vnl_matrix_fixed_ref_const<T,nrows,ncols>::div( const T* a, const T* b, T* r )
00085 {
00086   unsigned int count = nrows*ncols;
00087   while ( count-- )
00088     *(r++) = *(a++) / *(b++);
00089 }
00090 
00091 template<class T, unsigned nrows, unsigned ncols>
00092 void
00093 vnl_matrix_fixed_ref_const<T,nrows,ncols>::div( const T* a, T b, T* r )
00094 {
00095   unsigned int count = nrows*ncols;
00096   while ( count-- )
00097     *(r++) = *(a++) / b;
00098 }
00099 
00100 template<class T, unsigned nrows, unsigned ncols>
00101 bool
00102 vnl_matrix_fixed_ref_const<T,nrows,ncols>::equal( const T* a, const T* b )
00103 {
00104   unsigned int count = nrows*ncols;
00105   while ( count-- )
00106     if ( *(a++) != *(b++) )  return false;
00107   return true;
00108 }
00109 #endif // 0
00110 
00111 //------------------------------------------------------------
00112 
00113 
00114 template<class T, unsigned nrows, unsigned ncols>
00115 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00116 vnl_matrix_fixed_ref<T,nrows,ncols>::fill (T value) const
00117 {
00118   for (unsigned int i = 0; i < nrows; i++)
00119     for (unsigned int j = 0; j < ncols; j++)
00120       (*this)(i,j) = value;
00121   return *this;
00122 }
00123 
00124 
00125 template<class T, unsigned nrows, unsigned ncols>
00126 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00127 vnl_matrix_fixed_ref<T,nrows,ncols>::fill_diagonal(T value) const
00128 {
00129   for (unsigned int i = 0; i < nrows && i < ncols; i++)
00130     (*this)(i,i) = value;
00131   return *this;
00132 }
00133 
00134 
00135 template<class T, unsigned nrows, unsigned ncols>
00136 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00137 vnl_matrix_fixed_ref<T,nrows,ncols>::set_diagonal(vnl_vector<T> const& diag) const
00138 {
00139   assert(diag.size() >= nrows || diag.size() >= ncols);
00140   // The length of the diagonal of a non-square matrix is the minimum of
00141   // the matrix's width & height; that explains the "||" in the assert,
00142   // and the "&&" in the upper bound for the "for".
00143   for (unsigned int i = 0; i < nrows && i < ncols; ++i)
00144     (*this)(i,i) = diag[i];
00145   return *this;
00146 }
00147 
00148 
00149 template<class T, unsigned nrows, unsigned ncols>
00150 void
00151 vnl_matrix_fixed_ref_const<T,nrows,ncols>::print(vcl_ostream& os) const
00152 {
00153   for (unsigned int i = 0; i < nrows; i++)
00154   {
00155     os << (*this)(i,0);
00156     for (unsigned int j = 1; j < ncols; j++)
00157       os << ' ' << (*this)(i,j);
00158     os << '\n';
00159   }
00160 }
00161 
00162 
00163 template <class T, unsigned nrows, unsigned ncols>
00164 vnl_matrix_fixed<T,nrows,ncols>
00165 vnl_matrix_fixed_ref_const<T,nrows,ncols>::apply(T (*f)(T const&)) const
00166 {
00167   vnl_matrix_fixed<T,nrows,ncols> ret;
00168   vnl_c_vector<T>::apply(this->begin(), rows()*cols(), f, ret.data_block());
00169   return ret;
00170 }
00171 
00172 template <class T, unsigned nrows, unsigned ncols>
00173 vnl_matrix_fixed<T,nrows,ncols>
00174 vnl_matrix_fixed_ref_const<T,nrows,ncols>::apply(T (*f)(T)) const
00175 {
00176   vnl_matrix_fixed<T,nrows,ncols> ret;
00177   vnl_c_vector<T>::apply(this->begin(), rows()*cols(), f, ret.data_block());
00178   return ret;
00179 }
00180 
00181 ////--------------------------- Additions------------------------------------
00182 
00183 
00184 template<class T, unsigned nrows, unsigned ncols>
00185 vnl_matrix_fixed<T,ncols,nrows>
00186 vnl_matrix_fixed_ref_const<T,nrows,ncols>::transpose() const
00187 {
00188   vnl_matrix_fixed<T,ncols,nrows> result;
00189   for (unsigned int i = 0; i < cols(); i++)
00190     for (unsigned int j = 0; j < rows(); j++)
00191       result(i,j) = (*this)(j,i);
00192   return result;
00193 }
00194 
00195 template<class T, unsigned nrows, unsigned ncols>
00196 vnl_matrix_fixed<T,ncols,nrows>
00197 vnl_matrix_fixed_ref_const<T,nrows,ncols>::conjugate_transpose() const
00198 {
00199   vnl_matrix_fixed<T,ncols,nrows> result(transpose());
00200   vnl_c_vector<T>::conjugate(result.begin(),  // src
00201                              result.begin(),  // dst
00202                              result.size());  // size of block
00203   return result;
00204 }
00205 
00206 template<class T, unsigned nrows, unsigned ncols>
00207 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00208 vnl_matrix_fixed_ref<T,nrows,ncols>::update (vnl_matrix<T> const& m,
00209                                              unsigned top, unsigned left) const
00210 {
00211   const unsigned int bottom = top + m.rows();
00212   const unsigned int right = left + m.cols();
00213 #ifndef NDEBUG
00214   if (nrows < bottom || ncols < right)
00215     vnl_error_matrix_dimension ("update",
00216                                 bottom, right, m.rows(), m.cols());
00217 #endif
00218   for (unsigned int i = top; i < bottom; i++)
00219     for (unsigned int j = left; j < right; j++)
00220       (*this)(i,j) = m(i-top,j-left);
00221   return *this;
00222 }
00223 
00224 
00225 template<class T, unsigned nrows, unsigned ncols>
00226 vnl_matrix<T>
00227 vnl_matrix_fixed_ref_const<T,nrows,ncols>::extract (unsigned rowz, unsigned colz,
00228                                                     unsigned top, unsigned left) const
00229 {
00230 #ifndef NDEBUG
00231   unsigned int bottom = top + rowz;
00232   unsigned int right = left + colz;
00233   if ((nrows < bottom) || (ncols < right))
00234     vnl_error_matrix_dimension ("extract",
00235                                 nrows, ncols, bottom, right);
00236 #endif
00237   vnl_matrix<T> result(rowz, colz);
00238   for (unsigned int i = 0; i < rowz; i++)      // actual copy of all elements
00239     for (unsigned int j = 0; j < colz; j++)    // in submatrix
00240       result(i,j) = (*this)(top+i,left+j);
00241   return result;
00242 }
00243 
00244 
00245 template<class T, unsigned nrows, unsigned ncols>
00246 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00247 vnl_matrix_fixed_ref<T,nrows,ncols>::copy_in(T const *p) const
00248 {
00249   T* dp = this->data_block();
00250   unsigned int i = nrows * ncols;
00251   while (i--)
00252     *dp++ = *p++;
00253   return *this;
00254 }
00255 
00256 template<class T, unsigned nrows, unsigned ncols>
00257 void vnl_matrix_fixed_ref_const<T,nrows,ncols>::copy_out(T *p) const
00258 {
00259   T const* dp = this->data_block();
00260   unsigned int i = nrows*ncols;
00261   while (i--)
00262     *p++ = *dp++;
00263 }
00264 
00265 template<class T, unsigned nrows, unsigned ncols>
00266 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00267 vnl_matrix_fixed_ref<T,nrows,ncols>::set_identity() const
00268 {
00269   // Two simple loops are generally better than having a branch inside
00270   // the loop. Probably worth the O(n) extra writes.
00271   for (unsigned int i = 0; i < nrows; ++i)
00272     for (unsigned int j = 0; j < ncols; ++j)
00273         (*this)(i,j) = T(0);
00274   for (unsigned int i = 0; i < nrows && i < ncols; ++i)
00275     (*this)(i,i) = T(1);
00276   return *this;
00277 }
00278 
00279 //: Make each row of the matrix have unit norm.
00280 // All-zero rows are ignored.
00281 template<class T, unsigned nrows, unsigned ncols>
00282 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00283 vnl_matrix_fixed_ref<T,nrows,ncols>::normalize_rows() const
00284 {
00285   typedef typename vnl_numeric_traits<T>::abs_t Abs_t;
00286   for (unsigned int i = 0; i < nrows; i++)
00287   {
00288     Abs_t norm(0); // double will not do for all types.
00289     for (unsigned int j = 0; j < ncols; j++)
00290       norm += vnl_math_squared_magnitude( (*this)(i,j) );
00291 
00292     if (norm != 0)
00293     {
00294       typedef typename vnl_numeric_traits<Abs_t>::real_t real_t;
00295       real_t scale = real_t(1)/vcl_sqrt((real_t)norm);
00296       for (unsigned int j = 0; j < ncols; j++)
00297       {
00298         // FIXME need correct rounding here
00299         // There is e.g. no *standard* operator*=(complex<float>, double), hence the T() cast.
00300         (*this)(i,j) *= (T)(scale);
00301       }
00302     }
00303   }
00304   return *this;
00305 }
00306 
00307 template<class T, unsigned nrows, unsigned ncols>
00308 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00309 vnl_matrix_fixed_ref<T,nrows,ncols>::normalize_columns() const
00310 {
00311   typedef typename vnl_numeric_traits<T>::abs_t Abs_t;
00312   for (unsigned int j = 0; j < ncols; j++) {  // For each column in the Matrix
00313     Abs_t norm(0); // double will not do for all types.
00314     for (unsigned int i = 0; i < nrows; i++)
00315       norm += vnl_math_squared_magnitude( (*this)(i,j) );
00316 
00317     if (norm != 0)
00318     {
00319       typedef typename vnl_numeric_traits<Abs_t>::real_t real_t;
00320       real_t scale = real_t(1)/vcl_sqrt((real_t)norm);
00321       for (unsigned int i = 0; i < nrows; i++)
00322       {
00323         // FIXME need correct rounding here
00324         // There is e.g. no *standard* operator*=(complex<float>, double), hence the T() cast.
00325         (*this)(i,j) *= (T)(scale);
00326       }
00327     }
00328   }
00329   return *this;
00330 }
00331 
00332 template<class T, unsigned nrows, unsigned ncols>
00333 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00334 vnl_matrix_fixed_ref<T,nrows,ncols>::scale_row(unsigned row_index, T value) const
00335 {
00336 #ifndef NDEBUG
00337   if (row_index >= nrows)
00338     vnl_error_matrix_row_index("scale_row", row_index);
00339 #endif
00340   for (unsigned int j = 0; j < ncols; j++)
00341     (*this)(row_index,j) *= value;
00342   return *this;
00343 }
00344 
00345 template<class T, unsigned nrows, unsigned ncols>
00346 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00347 vnl_matrix_fixed_ref<T,nrows,ncols>::scale_column(unsigned column_index, T value) const
00348 {
00349 #ifndef NDEBUG
00350   if (column_index >= ncols)
00351     vnl_error_matrix_col_index("scale_column", column_index);
00352 #endif
00353   for (unsigned int j = 0; j < nrows; j++)
00354     (*this)(j,column_index) *= value;
00355   return *this;
00356 }
00357 
00358 //: Returns a copy of n rows, starting from "row"
00359 template<class T, unsigned nrows, unsigned ncols>
00360 vnl_matrix<T>
00361 vnl_matrix_fixed_ref_const<T,nrows,ncols>::get_n_rows (unsigned row, unsigned n) const
00362 {
00363 #ifndef NDEBUG
00364   if (row + n > nrows)
00365     vnl_error_matrix_row_index ("get_n_rows", row);
00366 #endif
00367 
00368   // Extract data rowwise.
00369   return vnl_matrix<T>((*this)[row], n, ncols);
00370 }
00371 
00372 template<class T, unsigned nrows, unsigned ncols>
00373 vnl_matrix<T>
00374 vnl_matrix_fixed_ref_const<T,nrows,ncols>::get_n_columns (unsigned column, unsigned n) const
00375 {
00376 #ifndef NDEBUG
00377   if (column + n > ncols)
00378     vnl_error_matrix_col_index ("get_n_columns", column);
00379 #endif
00380 
00381   vnl_matrix<T> result(nrows, n);
00382   for (unsigned int c = 0; c < n; ++c)
00383     for (unsigned int r = 0; r < nrows; ++r)
00384       result(r, c) = (*this)(r,column + c);
00385   return result;
00386 }
00387 
00388 #if 0 // commented out
00389 
00390 //: Create a vector out of row[row_index].
00391 template<class T, unsigned nrows, unsigned ncols>
00392 vnl_vector<T> vnl_matrix_fixed_ref_const<T,nrows,ncols>::get_row(unsigned row_index) const
00393 {
00394 #if ERROR_CHECKING
00395   if (row_index >= nrows)
00396     vnl_error_matrix_row_index ("get_row", row_index);
00397 #endif
00398 
00399   vnl_vector<T> v(ncols);
00400   for (unsigned int j = 0; j < ncols; j++)    // For each element in row
00401     v[j] = (*this)(row_index,j);
00402   return v;
00403 }
00404 
00405 //: Create a vector out of column[column_index].
00406 template<class T, unsigned nrows, unsigned ncols>
00407 vnl_vector<T> vnl_matrix_fixed_ref_const<T,nrows,ncols>::get_column(unsigned column_index) const
00408 {
00409 #if ERROR_CHECKING
00410   if (column_index >= ncols)
00411     vnl_error_matrix_col_index ("get_column", column_index);
00412 #endif
00413 
00414   vnl_vector<T> v(nrows);
00415   for (unsigned int j = 0; j < nrows; j++)
00416     v[j] = (*this)(j,column_index);
00417   return v;
00418 }
00419 
00420 #endif // 0
00421 
00422 //: Return a vector with the content of the (main) diagonal
00423 template<class T, unsigned nrows, unsigned ncols>
00424 vnl_vector<T> vnl_matrix_fixed_ref_const<T,nrows,ncols>::get_diagonal() const
00425 {
00426   vnl_vector<T> v(nrows < ncols ? nrows : ncols);
00427   for (unsigned int j = 0; j < nrows && j < ncols; ++j)
00428     v[j] = (*this)(j,j);
00429   return v;
00430 }
00431 
00432 //--------------------------------------------------------------------------------
00433 
00434 template<class T, unsigned nrows, unsigned ncols>
00435 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00436 vnl_matrix_fixed_ref<T,nrows,ncols>::set_row(unsigned row_index, T const *v) const
00437 {
00438   for (unsigned int j = 0; j < ncols; j++)
00439     (*this)(row_index,j) = v[j];
00440   return *this;
00441 }
00442 
00443 template<class T, unsigned nrows, unsigned ncols>
00444 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00445 vnl_matrix_fixed_ref<T,nrows,ncols>::set_row(unsigned row_index, vnl_vector_fixed<T,ncols> const &v) const
00446 {
00447   set_row(row_index,v.data_block());
00448   return *this;
00449 }
00450 
00451 template<class T, unsigned nrows, unsigned ncols>
00452 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00453 vnl_matrix_fixed_ref<T,nrows,ncols>::set_row(unsigned row_index, vnl_vector<T> const &v) const
00454 {
00455   set_row(row_index,v.data_block());
00456   return *this;
00457 }
00458 
00459 template<class T, unsigned nrows, unsigned ncols>
00460 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00461 vnl_matrix_fixed_ref<T,nrows,ncols>::set_row(unsigned row_index, T v) const
00462 {
00463   for (unsigned int j = 0; j < ncols; j++)
00464     (*this)(row_index,j) = v;
00465   return *this;
00466 }
00467 
00468 //--------------------------------------------------------------------------------
00469 
00470 template<class T, unsigned nrows, unsigned ncols>
00471 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00472 vnl_matrix_fixed_ref<T,nrows,ncols>::set_column(unsigned column_index, T const *v) const
00473 {
00474   for (unsigned int i = 0; i < nrows; i++)
00475     (*this)(i,column_index) = v[i];
00476   return *this;
00477 }
00478 
00479 template<class T, unsigned nrows, unsigned ncols>
00480 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00481 vnl_matrix_fixed_ref<T,nrows,ncols>::set_column(unsigned column_index, vnl_vector_fixed<T,nrows> const &v) const
00482 {
00483   set_column(column_index,v.data_block());
00484   return *this;
00485 }
00486 
00487 template<class T, unsigned nrows, unsigned ncols>
00488 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00489 vnl_matrix_fixed_ref<T,nrows,ncols>::set_column(unsigned column_index, vnl_vector<T> const &v) const
00490 {
00491   set_column(column_index,v.data_block());
00492   return *this;
00493 }
00494 
00495 template<class T, unsigned nrows, unsigned ncols>
00496 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00497 vnl_matrix_fixed_ref<T,nrows,ncols>::set_column(unsigned column_index, T v) const
00498 {
00499   for (unsigned int j = 0; j < nrows; j++)
00500     (*this)(j,column_index) = v;
00501   return *this;
00502 }
00503 
00504 
00505 template<class T, unsigned nrows, unsigned ncols>
00506 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00507 vnl_matrix_fixed_ref<T,nrows,ncols>::set_columns(unsigned starting_column, vnl_matrix<T> const& m) const
00508 {
00509 #ifndef NDEBUG
00510   if (nrows != m.rows() ||
00511       ncols < m.cols() + starting_column)           // Size match?
00512     vnl_error_matrix_dimension ("set_columns",
00513                                 nrows, ncols,
00514                                 m.rows(), m.cols());
00515 #endif
00516 
00517   for (unsigned int j = 0; j < m.cols(); ++j)
00518     for (unsigned int i = 0; i < nrows; i++)
00519       (*this)(i,starting_column + j) = m(i,j);
00520   return *this;
00521 }
00522 
00523 
00524 template <class T, unsigned nrows, unsigned ncols>
00525 bool
00526 vnl_matrix_fixed_ref_const<T,nrows,ncols>::is_identity() const
00527 {
00528   T const zero(0);
00529   T const one(1);
00530   for (unsigned int i = 0; i < nrows; ++i)
00531     for (unsigned int j = 0; j < ncols; ++j)
00532     {
00533       T xm = (*this)(i,j);
00534       if ( !((i == j) ? (xm == one) : (xm == zero)) )
00535         return false;
00536     }
00537   return true;
00538 }
00539 
00540 //: Return true if maximum absolute deviation of M from identity is <= tol.
00541 template <class T, unsigned nrows, unsigned ncols>
00542 bool
00543 vnl_matrix_fixed_ref_const<T,nrows,ncols>::is_identity(double tol) const
00544 {
00545   T one(1);
00546   for (unsigned int i = 0; i < nrows; ++i)
00547     for (unsigned int j = 0; j < ncols; ++j)
00548     {
00549       T xm = (*this)(i,j);
00550       abs_t absdev = (i == j) ? vnl_math_abs(xm - one) : vnl_math_abs(xm);
00551       if (absdev > tol)
00552         return false;
00553     }
00554   return true;
00555 }
00556 
00557 template <class T, unsigned nrows, unsigned ncols>
00558 bool
00559 vnl_matrix_fixed_ref_const<T,nrows,ncols>::is_zero() const
00560 {
00561   T const zero(0);
00562   for (unsigned int i = 0; i < nrows; ++i)
00563     for (unsigned int j = 0; j < ncols; ++j)
00564       if ( !( (*this)(i, j) == zero) )
00565         return false;
00566 
00567   return true;
00568 }
00569 
00570 template <class T, unsigned nrows, unsigned ncols>
00571 bool
00572 vnl_matrix_fixed_ref_const<T,nrows,ncols>::is_zero(double tol) const
00573 {
00574   for (unsigned int i = 0; i < nrows; ++i)
00575     for (unsigned int j = 0; j < ncols; ++j)
00576       if (vnl_math_abs((*this)(i,j)) > tol)
00577         return false;
00578 
00579   return true;
00580 }
00581 
00582 template <class T, unsigned nrows, unsigned ncols>
00583 bool
00584 vnl_matrix_fixed_ref_const<T,nrows,ncols>::has_nans() const
00585 {
00586   for (unsigned int i = 0; i < nrows; ++i)
00587     for (unsigned int j = 0; j < ncols; ++j)
00588       if (vnl_math_isnan((*this)(i,j)))
00589         return true;
00590 
00591   return false;
00592 }
00593 
00594 template <class T, unsigned nrows, unsigned ncols>
00595 bool
00596 vnl_matrix_fixed_ref_const<T,nrows,ncols>::is_finite() const
00597 {
00598   for (unsigned int i = 0; i < nrows; ++i)
00599     for (unsigned int j = 0; j < ncols; ++j)
00600       if (!vnl_math_isfinite((*this)(i,j)))
00601         return false;
00602 
00603   return true;
00604 }
00605 
00606 //: Abort if any element of M is inf or nan
00607 template <class T, unsigned nrows, unsigned ncols>
00608 void
00609 vnl_matrix_fixed_ref_const<T,nrows,ncols>::assert_finite_internal() const
00610 {
00611   if (is_finite())
00612     return;
00613 
00614   vcl_cerr << "\n\n" << __FILE__ " : " << __LINE__ << ": matrix has non-finite elements\n";
00615 
00616   if (rows() <= 20 && cols() <= 20)
00617     vcl_cerr << __FILE__ ": here it is:\n" << *this << '\n';
00618   else
00619   {
00620     vcl_cerr << __FILE__ ": it is quite big (" << rows() << 'x' << cols() << ")\n"
00621              << __FILE__ ": in the following picture '-' means finite and '*' means non-finite:\n";
00622 
00623     for (unsigned int i=0; i<rows(); ++i)
00624     {
00625       for (unsigned int j=0; j<cols(); ++j)
00626         vcl_cerr << char(vnl_math_isfinite((*this)(i, j)) ? '-' : '*');
00627       vcl_cerr << '\n';
00628     }
00629   }
00630   vcl_cerr << __FILE__ ": calling abort()\n";
00631   vcl_abort();
00632 }
00633 
00634 //: Abort unless M has the given size.
00635 template <class T, unsigned nrows, unsigned ncols>
00636 void
00637 vnl_matrix_fixed_ref_const<T,nrows,ncols>::assert_size_internal(unsigned rs,unsigned cs) const
00638 {
00639   if (nrows!=rs || ncols!=cs)
00640   {
00641     vcl_cerr << __FILE__ ": size is " << nrows << 'x' << ncols
00642              << ". should be " << rs << 'x' << cs << vcl_endl;
00643     vcl_abort();
00644   }
00645 }
00646 
00647 template <class T, unsigned nrows, unsigned ncols>
00648 bool
00649 vnl_matrix_fixed_ref<T,nrows,ncols>::read_ascii(vcl_istream& s) const
00650 {
00651   if (!s.good())
00652   {
00653     vcl_cerr << __FILE__ ": vnl_matrix_fixed_ref_const<T,nrows,ncols>::read_ascii: Called with bad stream\n";
00654     return false;
00655   }
00656 
00657   for (unsigned int i = 0; i < nrows; ++i)
00658     for (unsigned int j = 0; j < ncols; ++j)
00659       s >> (*this)(i,j);
00660 
00661   return s.good() || s.eof();
00662 }
00663 
00664 
00665 template <class T, unsigned nrows, unsigned ncols>
00666 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00667 vnl_matrix_fixed_ref<T,nrows,ncols>::flipud() const
00668 {
00669   for (unsigned int r1 = 0; 2*r1+1 < nrows; ++r1)
00670   {
00671     unsigned int r2 = nrows - 1 - r1;
00672     for (unsigned int c = 0; c < ncols; ++c)
00673     {
00674       T tmp = (*this)(r1, c);
00675       (*this)(r1, c) = (*this)(r2, c);
00676       (*this)(r2, c) = tmp;
00677     }
00678   }
00679   return *this;
00680 }
00681 
00682 
00683 template <class T, unsigned nrows, unsigned ncols>
00684 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00685 vnl_matrix_fixed_ref<T,nrows,ncols>::fliplr() const
00686 {
00687   for (unsigned int c1 = 0; 2*c1+1 < ncols; ++c1)
00688   {
00689     unsigned int c2 = ncols - 1 - c1;
00690     for (unsigned int r = 0; r < nrows; ++r)
00691     {
00692       T tmp = (*this)(r, c1);
00693       (*this)(r, c1) = (*this)(r, c2);
00694       (*this)(r, c2) = tmp;
00695     }
00696   }
00697   return *this;
00698 }
00699 
00700 template <class T, unsigned nrows, unsigned ncols>
00701 typename vnl_matrix_fixed_ref_const<T,nrows,ncols>::abs_t
00702 vnl_matrix_fixed_ref_const<T,nrows,ncols>::operator_one_norm() const
00703 {
00704   abs_t m(0);
00705   for (unsigned int j=0; j<ncols; ++j)
00706   {
00707     abs_t t(0);
00708     for (unsigned int i=0; i<nrows; ++i)
00709       t += vnl_math_abs( (*this)(i,j) );
00710     if (t > m)
00711       m = t;
00712   }
00713   return m;
00714 }
00715 
00716 template <class T, unsigned nrows, unsigned ncols>
00717 typename vnl_matrix_fixed_ref_const<T,nrows,ncols>::abs_t
00718 vnl_matrix_fixed_ref_const<T,nrows,ncols>::operator_inf_norm() const
00719 {
00720   abs_t m(0);
00721   for (unsigned int i=0; i<nrows; ++i)
00722   {
00723     abs_t t(0);
00724     for (unsigned int j=0; j<ncols; ++j)
00725       t += vnl_math_abs( (*this)(i,j) );
00726     if (t > m)
00727       m = t;
00728   }
00729   return m;
00730 }
00731 
00732 //: Transpose square matrix M in place.
00733 template <class T, unsigned nrows, unsigned ncols>
00734 vnl_matrix_fixed_ref<T,nrows,ncols> const&
00735 vnl_matrix_fixed_ref<T,nrows,ncols>::inplace_transpose() const
00736 {
00737   assert(nrows==ncols); // cannot inplace_transpose non-square fixed size matrix
00738   for (unsigned i = 0; i < nrows; ++i)
00739   for (unsigned j = i+1; j < ncols; ++j)
00740   {
00741     T t = (*this)(i,j);
00742     (*this)(i,j) = (*this)(j,i);
00743     (*this)(j,i) = t;
00744   }
00745   return *this;
00746 }
00747 
00748 
00749 #define VNL_MATRIX_FIXED_REF_INSTANTIATE(T,m,n) \
00750 template class vnl_matrix_fixed_ref_const<T, m, n >; \
00751 template class vnl_matrix_fixed_ref<T, m, n >
00752 
00753 #endif // vnl_matrix_fixed_ref_txx_