Basic Image AlgorithmS Library 2.8.0

fmat.h

00001 /*
00002 This file is distributed as part of the BIAS library (Basic ImageAlgorithmS)
00003 but it has not been developed by the authors of BIAS.
00004 
00005 For copyright, author and license information see below.
00006 */
00007 
00008 
00009 /*
00010 *
00011 * Template Numerical Toolkit (TNT): Linear Algebra Module
00012 *
00013 * Mathematical and Computational Sciences Division
00014 * National Institute of Technology,
00015 * Gaithersburg, MD USA
00016 *
00017 *
00018 * This software was developed at the National Institute of Standards and
00019 * Technology (NIST) by employees of the Federal Government in the course
00020 * of their official duties. Pursuant to title 17 Section 105 of the
00021 * United States Code, this software is not subject to copyright protection
00022 * and is in the public domain.  The Template Numerical Toolkit (TNT) is
00023 * an experimental system.  NIST assumes no responsibility whatsoever for
00024 * its use by other parties, and makes no guarantees, expressed or implied,
00025 * about its quality, reliability, or any other characteristic.
00026 *
00027 * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
00028 * see http://math.nist.gov/tnt for latest updates.
00029 *
00030 */
00031 
00032 
00033 
00034 // Fortran-compatible matrix: column oriented, 1-based (i,j) indexing
00035 
00036 #ifndef FMAT_H
00037 #define FMAT_H
00038 
00039 #include "subscript.h"
00040 #include "vec.h"
00041 #include <cstdlib>
00042 #include <cassert>
00043 #include <iostream>
00044 #include <sstream>
00045 #ifdef TNT_USE_REGIONS
00046 #include "region2d.h"
00047 #endif
00048 
00049 // simple 1-based, column oriented Matrix class
00050 
00051 namespace TNT
00052 {
00053 
00054   template <class T>
00055   class Fortran_Matrix 
00056   {
00057 
00058 
00059   public:
00060 
00061     typedef         T   value_type;
00062     typedef         T   element_type;
00063     typedef         T*  pointer;
00064     typedef         T*  iterator;
00065     typedef         T&  reference;
00066     typedef const   T*  const_iterator;
00067     typedef const   T&  const_reference;
00068 
00069     Subscript lbound() const { return 1;}
00070 
00071   protected:
00072     T* v_;                  // these are adjusted to simulate 1-offset
00073     Subscript m_;
00074     Subscript n_;
00075     T** col_;           // these are adjusted to simulate 1-offset
00076 
00077     // internal helper function to create the array
00078     // of row pointers
00079 
00080     void initialize(Subscript M, Subscript N)
00081     {
00082       // adjust col_[] pointers so that they are 1-offset:
00083       //   col_[j][i] is really col_[j-1][i-1];
00084       //
00085       // v_[] is the internal contiguous array, it is still 0-offset
00086       //
00087       v_ = new T[M*N];
00088       col_ = new T*[N];
00089 
00090       BIASASSERT(v_  != NULL);
00091       BIASASSERT(col_ != NULL);
00092 
00093 
00094       m_ = M;
00095       n_ = N;
00096       T* p = v_ - 1;              
00097       for (Subscript i=0; i<N; i++)
00098       {
00099         col_[i] = p;
00100         p += M ;
00101 
00102       }
00103       col_ --; 
00104     }
00105 
00106     void copy(const T*  v)
00107     {
00108       Subscript N = m_ * n_;
00109       Subscript i;
00110 
00111 #ifdef TNT_UNROLL_LOOPS
00112       Subscript Nmod4 = N & 3;
00113       Subscript N4 = N - Nmod4;
00114 
00115       for (i=0; i<N4; i+=4)
00116       {
00117         v_[i] = v[i];
00118         v_[i+1] = v[i+1];
00119         v_[i+2] = v[i+2];
00120         v_[i+3] = v[i+3];
00121       }
00122 
00123       for (i=N4; i< N; i++)
00124         v_[i] = v[i];
00125 #else
00126 
00127       for (i=0; i< N; i++)
00128         v_[i] = v[i];
00129 #endif      
00130     }
00131 
00132     void set(const T& val)
00133     {
00134       Subscript N = m_ * n_;
00135       Subscript i;
00136 
00137 #ifdef TNT_UNROLL_LOOPS
00138       Subscript Nmod4 = N & 3;
00139       Subscript N4 = N - Nmod4;
00140 
00141       for (i=0; i<N4; i+=4)
00142       {
00143         v_[i] = val;
00144         v_[i+1] = val;
00145         v_[i+2] = val;
00146         v_[i+3] = val; 
00147       }
00148 
00149       for (i=N4; i< N; i++)
00150         v_[i] = val;
00151 #else
00152 
00153       for (i=0; i< N; i++)
00154         v_[i] = val;
00155 
00156 #endif      
00157     }
00158 
00159 
00160 
00161     void destroy()
00162     {     
00163       /* do nothing, if no memory has been previously allocated */
00164       if (v_ == NULL) return ;
00165 
00166       /* if we are here, then matrix was previously allocated */
00167       delete [] (v_); v_ = NULL;    
00168       col_ ++;                // changed back to 0-offset
00169       delete [] (col_); col_ = NULL;
00170     }
00171 
00172 
00173   public:
00174 
00175     T* begin() { return v_; }
00176     const T* begin() const { return v_;}
00177 
00178     T* end() { return v_ + m_*n_; }
00179     const T* end() const { return v_ + m_*n_; }
00180 
00181 
00182     // constructors
00183 
00184     Fortran_Matrix() : v_(0), m_(0), n_(0), col_(0)  {};
00185     Fortran_Matrix(const Fortran_Matrix<T> &A)
00186     {
00187       initialize(A.m_, A.n_);
00188       copy(A.v_);
00189     }
00190 
00191     Fortran_Matrix(Subscript M, Subscript N, const T& value = T())
00192     {
00193       initialize(M,N);
00194       set(value);
00195     }
00196 
00197     Fortran_Matrix(Subscript M, Subscript N, const T* v)
00198     {
00199       initialize(M,N);
00200       copy(v);
00201     }
00202 
00203 
00204     Fortran_Matrix(Subscript M, Subscript N, char *s)
00205     {
00206       initialize(M,N);
00207       std::istringstream ins(s);
00208 
00209       Subscript i, j;
00210 
00211       for (i=1; i<=M; i++)
00212         for (j=1; j<=N; j++)
00213           ins >> (*this)(i,j);
00214     }
00215 
00216     // destructor
00217     ~Fortran_Matrix()
00218     {
00219       destroy();
00220     }
00221 
00222 
00223     // assignments
00224     //
00225     Fortran_Matrix<T>& operator=(const Fortran_Matrix<T> &A)
00226     {
00227       if (v_ == A.v_)
00228         return *this;
00229 
00230       if (m_ == A.m_  && n_ == A.n_)      // no need to re-alloc
00231         copy(A.v_);
00232 
00233       else
00234       {
00235         destroy();
00236         initialize(A.m_, A.n_);
00237         copy(A.v_);
00238       }
00239 
00240       return *this;
00241     }
00242 
00243     Fortran_Matrix<T>& operator=(const T& scalar)
00244     { 
00245       set(scalar); 
00246       return *this;
00247     }
00248 
00249 
00250     Subscript dim(Subscript d) const 
00251     {
00252 #ifdef TNT_BOUNDS_CHECK
00253       BIASASSERT( d >= 1);
00254       BIASASSERT( d <= 2);
00255 #endif
00256       return (d==1) ? m_ : ((d==2) ? n_ : 0); 
00257     }
00258 
00259     Subscript num_rows() const { return m_; }
00260     Subscript num_cols() const { return n_; }
00261 
00262     Fortran_Matrix<T>& newsize(Subscript M, Subscript N)
00263     {
00264       if (num_rows() == M && num_cols() == N)
00265         return *this;
00266 
00267       destroy();
00268       initialize(M,N);
00269 
00270       return *this;
00271     }
00272 
00273 
00274 
00275     // 1-based element access
00276     //
00277     inline reference operator()(Subscript i, Subscript j)
00278     { 
00279 #ifdef TNT_BOUNDS_CHECK
00280       BIASASSERT(1<=i);
00281       BIASASSERT(i <= m_) ;
00282       BIASASSERT(1<=j);
00283       BIASASSERT(j <= n_);
00284 #endif
00285       return col_[j][i]; 
00286     }
00287 
00288     inline const_reference operator() (Subscript i, Subscript j) const
00289     {
00290 #ifdef TNT_BOUNDS_CHECK
00291       BIASASSERT(1<=i);
00292       BIASASSERT(i <= m_) ;
00293       BIASASSERT(1<=j);
00294       BIASASSERT(j <= n_);
00295 #endif
00296       return col_[j][i]; 
00297     }
00298 
00299 
00300 #ifdef TNT_USE_REGIONS
00301 
00302     typedef Region2D<Fortran_Matrix<T> > Region;
00303     typedef const_Region2D< Fortran_Matrix<T> > const_Region;
00304 
00305     Region operator()(const Index1D &I, const Index1D &J)
00306     {
00307       return Region(*this, I,J);
00308     }
00309 
00310     const_Region operator()(const Index1D &I, const Index1D &J) const
00311     {
00312       return const_Region(*this, I,J);
00313     }
00314 
00315 #endif
00316 
00317 
00318   };
00319 
00320 
00321   /* ***************************  I/O  ********************************/
00322 
00323   template <class T>
00324     std::ostream& operator<<(std::ostream &s, const Fortran_Matrix<T> &A)
00325   {
00326     Subscript M=A.num_rows();
00327     Subscript N=A.num_cols();
00328 
00329     s << M << " " << N << "\n";
00330 
00331     for (Subscript i=1; i<=M; i++)
00332     {
00333       for (Subscript j=1; j<=N; j++)
00334       {
00335         s << A(i,j) << " ";
00336       }
00337       s << "\n";
00338     }
00339 
00340 
00341     return s;
00342   }
00343 
00344   template <class T>
00345     std::istream& operator>>(std::istream &s, Fortran_Matrix<T> &A)
00346   {
00347 
00348     Subscript M, N;
00349 
00350     s >> M >> N;
00351 
00352     if ( !(M == A.num_rows() && N == A.num_cols()))
00353     {
00354       A.newsize(M,N);
00355     }
00356 
00357 
00358     for (Subscript i=1; i<=M; i++)
00359       for (Subscript j=1; j<=N; j++)
00360       {
00361         s >>  A(i,j);
00362       }
00363 
00364 
00365       return s;
00366   }
00367 
00368   // *******************[ basic matrix algorithms ]***************************
00369 
00370 
00371   template <class T>
00372     Fortran_Matrix<T> operator+(const Fortran_Matrix<T> &A, 
00373     const Fortran_Matrix<T> &B)
00374   {
00375     Subscript M = A.num_rows();
00376     Subscript N = A.num_cols();
00377 
00378     BIASASSERT(M==B.num_rows());
00379     BIASASSERT(N==B.num_cols());
00380 
00381     Fortran_Matrix<T> tmp(M,N);
00382     Subscript i,j;
00383 
00384     for (i=1; i<=M; i++)
00385       for (j=1; j<=N; j++)
00386         tmp(i,j) = A(i,j) + B(i,j);
00387 
00388     return tmp;
00389   }
00390 
00391   template <class T>
00392     Fortran_Matrix<T> operator-(const Fortran_Matrix<T> &A, 
00393     const Fortran_Matrix<T> &B)
00394   {
00395     Subscript M = A.num_rows();
00396     Subscript N = A.num_cols();
00397 
00398     BIASASSERT(M==B.num_rows());
00399     BIASASSERT(N==B.num_cols());
00400 
00401     Fortran_Matrix<T> tmp(M,N);
00402     Subscript i,j;
00403 
00404     for (i=1; i<=M; i++)
00405       for (j=1; j<=N; j++)
00406         tmp(i,j) = A(i,j) - B(i,j);
00407 
00408     return tmp;
00409   }
00410 
00411   // element-wise multiplication  (use matmult() below for matrix
00412   // multiplication in the linear algebra sense.)
00413   //
00414   //
00415   template <class T>
00416     Fortran_Matrix<T> mult_element(const Fortran_Matrix<T> &A, 
00417     const Fortran_Matrix<T> &B)
00418   {
00419     Subscript M = A.num_rows();
00420     Subscript N = A.num_cols();
00421 
00422     BIASASSERT(M==B.num_rows());
00423     BIASASSERT(N==B.num_cols());
00424 
00425     Fortran_Matrix<T> tmp(M,N);
00426     Subscript i,j;
00427 
00428     for (i=1; i<=M; i++)
00429       for (j=1; j<=N; j++)
00430         tmp(i,j) = A(i,j) * B(i,j);
00431 
00432     return tmp;
00433   }
00434 
00435 
00436   template <class T>
00437     Fortran_Matrix<T> transpose(const Fortran_Matrix<T> &A)
00438   {
00439     Subscript M = A.num_rows();
00440     Subscript N = A.num_cols();
00441 
00442     Fortran_Matrix<T> S(N,M);
00443     Subscript i, j;
00444 
00445     for (i=1; i<=M; i++)
00446       for (j=1; j<=N; j++)
00447         S(j,i) = A(i,j);
00448 
00449     return S;
00450   }
00451 
00452 
00453 
00454   template <class T>
00455     inline Fortran_Matrix<T> matmult(const Fortran_Matrix<T>  &A, 
00456     const Fortran_Matrix<T> &B)
00457   {
00458 
00459 #ifdef TNT_BOUNDS_CHECK
00460     BIASASSERT(A.num_cols() == B.num_rows());
00461 #endif
00462 
00463     Subscript M = A.num_rows();
00464     Subscript N = A.num_cols();
00465     Subscript K = B.num_cols();
00466 
00467     Fortran_Matrix<T> tmp(M,K);
00468     T sum;
00469 
00470     for (Subscript i=1; i<=M; i++)
00471       for (Subscript k=1; k<=K; k++)
00472       {
00473         sum = 0;
00474         for (Subscript j=1; j<=N; j++)
00475           sum = sum +  A(i,j) * B(j,k);
00476 
00477         tmp(i,k) = sum; 
00478       }
00479 
00480       return tmp;
00481   }
00482 
00483   template <class T>
00484     inline Fortran_Matrix<T> operator*(const Fortran_Matrix<T> &A, 
00485     const Fortran_Matrix<T> &B)
00486   {
00487     return matmult(A,B);
00488   }
00489 
00490   template <class T>
00491     inline int matmult(Fortran_Matrix<T>& C, const Fortran_Matrix<T>  &A, 
00492     const Fortran_Matrix<T> &B)
00493   {
00494 
00495     BIASASSERT(A.num_cols() == B.num_rows());
00496 
00497     Subscript M = A.num_rows();
00498     Subscript N = A.num_cols();
00499     Subscript K = B.num_cols();
00500 
00501     C.newsize(M,K);         // adjust shape of C, if necessary
00502 
00503 
00504     T sum; 
00505 
00506     const T* row_i;
00507     const T* col_k;
00508 
00509     for (Subscript i=1; i<=M; i++)
00510     {
00511       for (Subscript k=1; k<=K; k++)
00512       {
00513         row_i = &A(i,1);
00514         col_k = &B(1,k);
00515         sum = 0;
00516         for (Subscript j=1; j<=N; j++)
00517         {
00518           sum +=  *row_i * *col_k;
00519           row_i += M;
00520           col_k ++;
00521         }
00522 
00523         C(i,k) = sum; 
00524       }
00525 
00526     }
00527 
00528     return 0;
00529   }
00530 
00531 
00532   template <class T>
00533     Vector<T> matmult(const Fortran_Matrix<T>  &A, const Vector<T> &x)
00534   {
00535 
00536 #ifdef TNT_BOUNDS_CHECK
00537     BIASASSERT(A.num_cols() == x.dim());
00538 #endif
00539 
00540     Subscript M = A.num_rows();
00541     Subscript N = A.num_cols();
00542 
00543     Vector<T> tmp(M);
00544     T sum;
00545 
00546     for (Subscript i=1; i<=M; i++)
00547     {
00548       sum = 0;
00549       for (Subscript j=1; j<=N; j++)
00550         sum = sum +  A(i,j) * x(j);
00551 
00552       tmp(i) = sum; 
00553     }
00554 
00555     return tmp;
00556   }
00557 
00558   template <class T>
00559     inline Vector<T> operator*(const Fortran_Matrix<T>  &A, const Vector<T> &x)
00560   {
00561     return matmult(A,x);
00562   }
00563 
00564   template <class T>
00565     inline Fortran_Matrix<T> operator*(const Fortran_Matrix<T>  &A, const T &x)
00566   {
00567     Subscript M = A.num_rows();
00568     Subscript N = A.num_cols();
00569 
00570     Subscript MN = M*N; 
00571 
00572     Fortran_Matrix<T> res(M,N);
00573     const T* a = A.begin();
00574     T* t = res.begin();
00575     T* tend = res.end();
00576 
00577     for (t=res.begin(); t < tend; t++, a++)
00578       *t = *a * x;
00579 
00580     return res;
00581   } 
00582 
00583 }  // namespace TNT
00584 #endif
00585 // FMAT_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends