Main Page | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members

smatrix.h

Go to the documentation of this file.
00001 // -*-c++-*-
00002 #ifndef IG_SMATRIX_H
00003 #define IG_SMATRIX_H
00004 
00005 // $Id: smatrix.h,v 1.29 2005/07/02 15:58:26 hkuiper Exp $
00006 
00007 // CwMtx matrix and vector math library
00008 // Copyright (C) 1999-2001  Harry Kuiper
00009 // Copyright (C) 2000  Will DeVore (template conversion)
00010 
00011 // This library is free software; you can redistribute it and/or
00012 // modify it under the terms of the GNU Lesser General Public
00013 // License as published by the Free Software Foundation; either
00014 // version 2 of the License, or (at your option) any later version.
00015 
00016 // This library is distributed in the hope that it will be useful,
00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 
00026 #ifndef IG_MATRIX_H
00027 #include "matrix.h"
00028 #endif
00029 
00030 namespace CwMtx
00031 {
00032   // prefix smat
00033   template < class T = double >
00034   class CWTSquareMatrix : public CWTMatrix<T>
00035   {
00036   public:
00037     CWTSquareMatrix(): CWTMatrix<T>() {};
00038     CWTSquareMatrix(unsigned crowInit): CWTMatrix<T>( crowInit, crowInit) {};
00039     CWTSquareMatrix(const CWTMatrix<T> &mat): CWTMatrix<T>(mat) {};
00040     CWTSquareMatrix(const CWTSquareMatrix &smat): CWTMatrix<T>(smat) {};
00041     CWTSquareMatrix(const CWTMatrix<T> &, unsigned, unsigned, unsigned);
00042     CWTSquareMatrix(const CWTSquareMatrix<T> &, unsigned, unsigned, unsigned);
00043 
00044     ~CWTSquareMatrix() {};
00045 
00046     void dimension(unsigned crowInit)
00047     {
00048       CWTMatrix<T>::dimension(crowInit, crowInit);
00049     }
00050 
00051     void mapInto(const CWTSquareMatrix &, unsigned, unsigned, unsigned);
00052 
00053     CWTSquareMatrix operator +(const CWTSquareMatrix &) const;
00054     CWTSquareMatrix operator -(const CWTSquareMatrix &) const;
00055     CWTSquareMatrix operator -() const;
00056     CWTSquareMatrix operator *(const T &) const;
00057     CWTSquareMatrix operator *(const CWTSquareMatrix &) const;
00058     CWTSquareMatrix operator /(const T &value) const
00059     {
00060       // NOTE: This used to work with g++ <= 3.2.
00061       // return (*this)*static_cast<const T &>(CWTUnity<T>()/value);
00062 
00063       T tTmp = CWTUnity<T>();
00064       tTmp /= value;
00065       return (*this)*tTmp;
00066     }
00067     CWTSquareMatrix operator /(const CWTSquareMatrix &) const;
00068 
00069     // not inherited
00070     CWTSquareMatrix & operator =(const CWTSquareMatrix &smat);
00071     CWTSquareMatrix & operator +=(const CWTSquareMatrix &smat);
00072     CWTSquareMatrix & operator -=(const CWTSquareMatrix &smat);
00073     CWTSquareMatrix & operator *=(const T &value);
00074     CWTSquareMatrix & operator *=(const CWTSquareMatrix &);
00075     CWTSquareMatrix & operator /=(const T &value);
00076     CWTSquareMatrix & operator /=(const CWTSquareMatrix &);
00077 
00078     // stores the adjoint of argument in this
00079     void storeAdjoint(const CWTSquareMatrix &);
00080     // stores the inverse of argument in this
00081     void storeInverse(const CWTSquareMatrix &);
00082     // makes this its own adjoint
00083     void makeAdjoint();
00084     // makes this its own inverse
00085     void makeInverse();
00086     // makes this a unity matrix
00087     void makeUnity();
00088   };
00089 
00090   template <class T, unsigned crow>
00091   class CWTSMat: public T
00092   {
00093   public:
00094     CWTSMat(): T(crow) {}
00095 
00096     T & operator =(const T &smtx) { return T::operator=(smtx); }
00097   };
00098 
00099   // Unity square matrix
00100   template <class T, unsigned crow>
00101   class CWTUnity< CWTSMat<CWTSquareMatrix<T>, crow> >:
00102     public CWTSMat<CWTSquareMatrix<T>, crow>
00103   {
00104   public:
00105     CWTUnity() { this->makeUnity(); }
00106   };
00107 
00108   // Zero matrix.
00109   template <class T, unsigned crow>
00110   class CWTZero< CWTSMat<CWTSquareMatrix<T>, crow> >:
00111     public CWTSMat<CWTSquareMatrix<T>, crow>
00112   {
00113   public:
00114     CWTZero() { fill(CWTZero<T>()); }
00115   };
00116 
00117   //
00118   // Constructors
00119   //
00120 
00121   template < class T >
00122   inline CWTSquareMatrix<T>::CWTSquareMatrix(const CWTMatrix<T> &mat,
00123                                              unsigned irowStart,
00124                                              unsigned icolStart,
00125                                              unsigned irowEnd)
00126     :
00127     CWTMatrix<T>(mat,
00128                  irowStart,
00129                  icolStart,
00130                  irowEnd,
00131                  icolStart + irowEnd - irowStart)
00132   {
00133   }
00134 
00135   template < class T >
00136   inline CWTSquareMatrix<T>::CWTSquareMatrix(const CWTSquareMatrix<T> &smat,
00137                                              unsigned irowStart,
00138                                              unsigned icolStart,
00139                                              unsigned irowEnd)
00140     :
00141     CWTMatrix<T>(smat,
00142                  irowStart,
00143                  icolStart,
00144                  irowEnd,
00145                  icolStart + irowEnd - irowStart)
00146   {
00147   }
00148 
00149   //
00150   // User Methods
00151   //
00152 
00153   template < class T >
00154   inline void CWTSquareMatrix<T>::mapInto(const CWTSquareMatrix<T> &smat,
00155                                           unsigned irowStart,
00156                                           unsigned icolStart,
00157                                           unsigned irowEnd)
00158   {
00159     CWTMatrix<T>::mapInto(smat,
00160                           irowStart,
00161                           icolStart,
00162                           irowEnd,
00163                           icolStart + irowEnd - irowStart);
00164   }
00165 
00166   // not inherited
00167   template < class T >
00168   inline CWTSquareMatrix<T> &
00169   CWTSquareMatrix<T>::operator =(const CWTSquareMatrix<T> &smat)
00170   {
00171     return static_cast<CWTSquareMatrix<T> &>(CWTMatrix<T>::operator=(smat));
00172   }
00173 
00174   template < class T >
00175   inline CWTSquareMatrix<T> &
00176   CWTSquareMatrix<T>::operator +=(const CWTSquareMatrix<T> &smat)
00177   {
00178     return static_cast<CWTSquareMatrix<T> &>(CWTMatrix<T>::operator+=(smat));
00179   }
00180 
00181   template < class T >
00182   inline CWTSquareMatrix<T> &
00183   CWTSquareMatrix<T>::operator -=(const CWTSquareMatrix &smat)
00184   {
00185     return static_cast<CWTSquareMatrix<T> &>(CWTMatrix<T>::operator-=(smat));
00186   }
00187 
00188   template < class T >
00189   inline CWTSquareMatrix<T> &
00190   CWTSquareMatrix<T>::operator *=(const T &value)
00191   {
00192     return static_cast<CWTSquareMatrix<T> &>(CWTMatrix<T>::operator*=(value));
00193   }
00194 
00195   template < class T >
00196   inline CWTSquareMatrix<T> &
00197   CWTSquareMatrix<T>::operator *=(const CWTSquareMatrix<T> &smat)
00198   {
00199     return (*this) = (*this)*smat;
00200   }
00201 
00202   template < class T >
00203   inline CWTSquareMatrix<T> &
00204   CWTSquareMatrix<T>::operator /=(const T &value)
00205   {
00206     // NOTE: This used to work with g++ <= 3.2.
00207     // return (*this) *= static_cast<const T &>(CWTUnity<T>()/value);
00208 
00209     T tTmp = CWTUnity<T>();
00210     tTmp /= value;
00211     return (*this) *= tTmp;
00212   }
00213 
00214   template < class T >
00215   inline CWTSquareMatrix<T> &
00216   CWTSquareMatrix<T>::operator /=(const CWTSquareMatrix<T> &smat)
00217   {
00218     return (*this) = (*this)/smat;
00219   }
00220 
00221   template < class T >
00222   CWTSquareMatrix<T>
00223   CWTSquareMatrix<T>::operator +(const CWTSquareMatrix<T> &smat) const
00224   {
00225     return CWTSquareMatrix(*this) += smat;
00226   }
00227 
00228   template < class T >
00229   CWTSquareMatrix<T>
00230   CWTSquareMatrix<T>::operator -(const CWTSquareMatrix<T> &smat) const
00231   {
00232     return CWTSquareMatrix(*this) -= smat;
00233   }
00234 
00235   template < class T >
00236   CWTSquareMatrix<T> CWTSquareMatrix<T>::operator -() const
00237   {
00238     // NOTE: This used to work with g++ <= 3.2.
00239     // return (*this)*static_cast<const T &>(CWTZero<T>() - CWTUnity<T>());
00240 
00241     T tTmp = CWTZero<T>();
00242     tTmp -= CWTUnity<T>();
00243     return (*this)*tTmp;
00244   }
00245 
00246   template < class T >
00247   CWTSquareMatrix<T> CWTSquareMatrix<T>::operator *(const T &value) const
00248   {
00249     return CWTSquareMatrix(*this) *= value;
00250   }
00251 
00252   template < class T >
00253   CWTSquareMatrix<T>
00254   CWTSquareMatrix<T>::operator *(const CWTSquareMatrix<T> &smat) const
00255   {
00256     CWTSquareMatrix smatResult(this->getRows());
00257     smatResult.storeProduct(*this, smat);
00258     return smatResult;
00259   }
00260 
00261   template < class T >
00262   CWTSquareMatrix<T>
00263   CWTSquareMatrix<T>::operator /(const CWTSquareMatrix<T> &smat) const
00264   {
00265     return (*this)*inv(smat);
00266   }
00267 
00268   // stores inverse of argument in this
00269   template < class T >
00270   void CWTSquareMatrix<T>::storeInverse(const CWTSquareMatrix<T> &smat)
00271   {
00272     // copy input matrix in this
00273     (*this) = smat;
00274     makeInverse();
00275   }
00276 
00277   // makes this a unity matrix
00278   template < class T >
00279   void CWTSquareMatrix<T>::makeUnity()
00280   {
00281     unsigned crow(this->getRows()), ccol(this->getCols());
00282 
00283     for (unsigned irow = 0; irow < crow; ++irow)
00284       {
00285         for (unsigned icol = 0; icol < ccol; ++icol)
00286           {
00287             if (irow == icol)
00288               {
00289                 (*this)[irow][icol] = CWTUnity<T>();
00290               }
00291             else
00292               {
00293                 (*this)[irow][icol] = CWTZero<T>();
00294               }
00295           }
00296       }
00297   }
00298 
00299   // makes this its own adjoint
00300   template < class T >
00301   void CWTSquareMatrix<T>::makeAdjoint()
00302   {
00303     // we need a copy of this
00304     CWTSquareMatrix smatCopy(*this);
00305     // for easier access to crows
00306     unsigned crowCopy = smatCopy.getRows();
00307     T elemTmp;
00308     // will eventually contain det(smatCopy)
00309     T elemDet = CWTUnity<T>();
00310 
00311     // make this a unity matrix
00312     makeUnity();
00313 
00314     // start row reduction
00315     for (unsigned irow = 0; irow < crowCopy; ++irow)
00316       {
00317         // if element [irow][irow] is zero, add a row with non-zero
00318         // element at column irow
00319         if (smatCopy[irow][irow] == CWTZero<T>())
00320           {
00321             for (unsigned irowTmp = irow; irowTmp < crowCopy; ++irowTmp)
00322               {
00323                 if (smatCopy[irowTmp][irow] != CWTZero<T>())
00324                   {
00325                     // has no effect on adj(smat)
00326                     smatCopy.addRowToRow(irowTmp, irow);
00327                     // repeat action on this
00328                     this->addRowToRow(irowTmp, irow);
00329                     break;
00330                   };
00331               };
00332           };
00333 
00334         elemTmp = smatCopy[irow][irow];
00335         T tTmp = CWTUnity<T>();
00336         tTmp /= elemTmp;
00337         smatCopy.multiplyRow(irow, tTmp);
00338         // repeat action on this
00339         multiplyRow(irow, tTmp);
00340         elemDet *= elemTmp;
00341 
00342         for (unsigned irowToAddTo = 0; irowToAddTo < crowCopy; ++irowToAddTo)
00343           {
00344             if (irowToAddTo != irow)
00345               {
00346                 elemTmp = -smatCopy[irowToAddTo][irow];
00347                 smatCopy.addRowToRow(irow, irowToAddTo, elemTmp);
00348                 // repeat action on this
00349                 addRowToRow(irow, irowToAddTo, elemTmp);
00350               };
00351           };
00352       };
00353 
00354     // this now contains its adjoint
00355     (*this) *= elemDet;
00356   }
00357 
00358   // stores adjoint of input matrix in this
00359   template < class T >
00360   void CWTSquareMatrix<T>::storeAdjoint(const CWTSquareMatrix<T> &smat)
00361   {
00362     // copy input matrix in this
00363     (*this) = smat;
00364     makeAdjoint();
00365   }
00366 
00367   // makes this its own inverse
00368   template < class T >
00369   void CWTSquareMatrix<T>::makeInverse()
00370   {
00371     // we need a copy of this
00372     CWTSquareMatrix smatCopy(*this);
00373     // for easier access to crows
00374     unsigned crowCopy = smatCopy.getRows();
00375     T elemTmp;
00376 
00377     // make this a unity matrix
00378     makeUnity();
00379 
00380     // start row reduction
00381     for (unsigned irow = 0; irow < crowCopy; ++irow)
00382       {
00383         // if element [irow][irow] is zero, add a row with non-zero
00384         // element at column irow
00385         if (smatCopy[irow][irow] == CWTZero<T>())
00386           {
00387             for (unsigned irowTmp = irow; irowTmp < crowCopy; ++irowTmp)
00388               {
00389                 // has no effect on inv(smat)
00390                 if (smatCopy[irowTmp][irow] != CWTZero<T>())
00391                   {
00392                     smatCopy.addRowToRow(irowTmp, irow);
00393                     // repeat action on this
00394                     this->addRowToRow(irowTmp, irow);
00395                     break;
00396                   };
00397               };
00398           };
00399 
00400         elemTmp = smatCopy[irow][irow];
00401 
00402         // NOTE: This used to work with g++ <= 3.2.
00403         // smatCopy.multiplyRow(irow,
00404         //                   static_cast<const T &>(CWTUnity<T>()/elemTmp));
00405         // multiplyRow(irow, static_cast<const T &>(CWTUnity<T>()/elemTmp));
00406 
00407         T tTmp = CWTUnity<T>();
00408         tTmp /= elemTmp;
00409         smatCopy.multiplyRow(irow, tTmp);
00410         multiplyRow(irow, tTmp);
00411 
00412         for (unsigned irowToAddTo = 0; irowToAddTo < crowCopy; ++irowToAddTo)
00413           {
00414             if (irowToAddTo != irow)
00415               {
00416                 elemTmp = -smatCopy[irowToAddTo][irow];
00417                 smatCopy.addRowToRow(irow, irowToAddTo, elemTmp);
00418                 // repeat action on this
00419                 addRowToRow(irow, irowToAddTo, elemTmp);
00420               };
00421           };
00422       };
00423 
00424     // this now contains its inverse
00425   }
00426 
00427   //
00428   // template functions designed work with the Square Matrix class.
00429   //
00430 
00431   // SCLR*CWTSquareMatrix
00432   template < class T >
00433   inline CWTSquareMatrix<T> operator *(const T &value,
00434                                        const CWTSquareMatrix<T> &smat)
00435   {
00436     return smat*value;
00437   }
00438 
00439   template < class T >
00440   CWTSquareMatrix<T> transpose(const CWTSquareMatrix<T> &smat)
00441   {
00442     CWTSquareMatrix<T> smatResult(smat.getRows());
00443     smatResult.storeTranspose(smat);
00444     return smatResult;
00445   }
00446 
00447   // returns the adjoint of a square matrix
00448   template < class T >
00449   CWTSquareMatrix<T> adj(const CWTSquareMatrix<T> &smat)
00450   {
00451     CWTSquareMatrix<T> smatResult(smat);   // copy input matrix
00452     smatResult.makeAdjoint();
00453     return smatResult;
00454   }
00455 
00456   // calculates the inverse of a square matrix
00457   template < class T >
00458   CWTSquareMatrix<T> inv(const CWTSquareMatrix<T> &smat)
00459   {
00460     // copy input matrix
00461     CWTSquareMatrix<T> smatResult(smat);
00462     smatResult.makeInverse();
00463     return smatResult;
00464   }
00465 
00466   // calculates the determinant of a square matrix
00467   template < class T >
00468   T det(const CWTSquareMatrix<T> &smat)
00469   {
00470     // a copy of the input matrix
00471     CWTSquareMatrix<T> smatCopy(smat);
00472     unsigned crowCopy = smatCopy.getRows();
00473 
00474     // start row reduction
00475     T elemTmp;
00476     // will eventually contain det(smatCopy)
00477     T elemDet = CWTUnity<T>();
00478 
00479     for (unsigned irow = 0; irow < crowCopy; ++irow)
00480       {
00481         // if element [irow][irow] is zero, add a row with non-zero
00482         // element at column irow
00483         if (smatCopy[irow][irow] == CWTZero<T>())
00484           {
00485             for (unsigned irowTmp = irow; irowTmp < crowCopy; ++irowTmp)
00486               {
00487                 // has no effect on inv(smat)
00488                 if (smatCopy[irowTmp][irow] != CWTZero<T>())
00489                   {
00490                     smatCopy.addRowToRow(irowTmp, irow);
00491                     break;
00492                   };
00493               };
00494           };
00495 
00496         elemTmp =  smatCopy[irow][irow];
00497         elemDet *= elemTmp;
00498 
00499         if (elemDet == CWTZero<T>())
00500           {
00501             // once elemDet is zero it will stay zero
00502             return elemDet;
00503           }
00504 
00505         smatCopy.multiplyRow(irow, CWTUnity<T>()/elemTmp);
00506 
00507         for (unsigned irowToAddTo = 0; irowToAddTo < crowCopy; ++irowToAddTo)
00508           {
00509             if (irowToAddTo != irow)
00510               {
00511                 smatCopy.addRowToRow(irow,
00512                                      irowToAddTo,
00513                                      -smatCopy[irowToAddTo][irow]);
00514               };
00515           };
00516       };
00517 
00518     return elemDet;
00519   }
00520 
00521   // trace
00522   template < class T >
00523   T tr(const CWTSquareMatrix<T> &smat)
00524   {
00525     T elemTmp = CWTZero<T>();
00526 
00527     for (unsigned c = 0; c < smat.getCols(); c++)
00528       {
00529         elemTmp += smat[c][c];
00530       }
00531 
00532     return elemTmp;
00533   }
00534 
00535 }
00536 
00537 #endif // IG_SMATRIX_H

Generated on Sun Jul 3 12:18:14 2005 for Matrix and vector library by  doxygen 1.4.2