LocalMatrix.h

Go to the documentation of this file.
00001 /*==============================================================================
00002 
00003                                     O  F  E  L  I
00004 
00005                            Object  Finite  Element  Library
00006 
00007   ==============================================================================
00008 
00009    Copyright (C) 1998 - 2008 Rachid Touzani
00010 
00011    This program is free software; you can redistribute it and/or modify it under
00012    the terms of the GNU General Public License as published by the Free
00013    Software Foundation; Version 2 of the License.
00014 
00015    This program is distributed in the hope that it will be useful, but WITHOUT
00016    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00017    FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
00018    details.
00019 
00020    You should have received a copy of the GNU General Public License
00021    along with this program; if not, write to the :
00022 
00023    Free Software Foundation
00024    Inc., 59 Temple Place - Suite 330
00025    Boston, MA  02111-1307, USA
00026 
00027   ==============================================================================
00028 
00029       Definition of Template Class LocalMatrix for dense small size matrices
00030 
00031   ==============================================================================*/
00032 
00033 
00034 #ifndef __LOCAL_MATRIX_H
00035 #define __LOCAL_MATRIX_H
00036 
00037 
00038 #include "OFELI_Config.h"
00039 #include "LocalVect.h"
00040 #include "SkMatrix.h"
00041 #include "SkSMatrix.h"
00042 #include "SpMatrix.h"
00043 
00044 namespace OFELI {
00045 
00066 template<class T_, size_t NC_> class LocalVect;
00067 template<class T_> class SkMatrix;
00068 template<class T_> class SkSMatrix;
00069 template<class T_> class SpMatrix;
00070 
00071 
00072 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00073 template<class T_, size_t NR_, size_t NC_>
00074 struct ErrorInLocalMatrix {
00075    void Message(const char *file, size_t line, int code, int p1=0, int p2=0, T_ p3=T_(0))
00076    {
00077      cerr << "\n*** Fatal Error in OFELI ***\n";
00078      cerr << "File : " << file << ", line : " << line << "\nIn LocalMatrix<>::";
00079      switch (code) {
00080 
00081         case 12:
00082            cerr << "Factor() : The " << p1 << "-th pivot = " << p3 << " is too small\n";
00083            break;
00084 
00085         case 11:
00086            cerr << "Factor() : Can't factorize a rectangle matrix\n";
00087            break;
00088 
00089         case 21:
00090            cerr << "Solve(b) : The " << p1 << "-th diagonal entry = " << p3 << " is too small\n";
00091            break;
00092 
00093         case 22:
00094            cerr << "Solve(b) : Can't solve with a rectangle matrix\n";
00095 
00096         case 31:
00097            cerr << "set() : This argument is valid for square matrices only\n";
00098      }
00099      cerr << "Program stops" << endl;
00100      exit(1);
00101    }
00102 };
00103 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00104 
00105 
00106 template<class T_, size_t NR_, size_t NC_> class LocalMatrix
00107 {
00108 
00109  public :
00110 
00113     LocalMatrix();
00114 
00116     LocalMatrix(const LocalMatrix<T_,NR_,NC_> &m);
00117 
00122     LocalMatrix(Element *el, const SpMatrix<T_> &a);
00123 
00128     LocalMatrix(Element *el, const SkMatrix<T_> &a);
00129 
00134     LocalMatrix(Element *el, const SkSMatrix<T_> &a);
00135 
00137     ~LocalMatrix();
00138 
00141     T_ & operator()(size_t i, size_t j) { return _a[(i-1)*NC_+j-1]; }
00142 
00145     T_ operator()(size_t i, size_t j) const { return _a[(i-1)*NC_+j-1]; }
00146 
00147 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00148     void set(int opt);
00149 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00150 
00156     void Localize(Element *el, const SpMatrix<T_> &a);
00157 
00163     void Localize(Element *el, const SkMatrix<T_> &a);
00164 
00170     void Localize(Element *el, const SkSMatrix<T_> &a);
00171 
00174     LocalMatrix<T_,NR_,NC_> & operator=(const LocalMatrix<T_,NR_,NC_> &m);
00175 
00178     LocalMatrix<T_,NR_,NC_> & operator=(const T_ &x);
00179 
00182     LocalMatrix<T_,NR_,NC_> & operator+=(const LocalMatrix<T_,NR_,NC_> &m);
00183 
00186     LocalMatrix<T_,NR_,NC_> & operator-=(const LocalMatrix<T_,NR_,NC_> &m);
00187 
00190     LocalVect<T_,NR_> operator*(LocalVect<T_,NC_> &x);
00191 
00194     LocalMatrix<T_,NR_,NC_> & operator+=(const T_ &x);
00195 
00198     LocalMatrix<T_,NR_,NC_> & operator-=(const T_ &x);
00199 
00202     LocalMatrix<T_,NR_,NC_> & operator*=(const T_ &x);
00203 
00206     LocalMatrix<T_,NR_,NC_> & operator/=(const T_ &x);
00207 
00212     void MultAdd(const LocalVect<T_,NC_> &x, LocalVect<T_,NR_> &y);
00213 
00219     void MultAddScal(const T_ &a, const LocalVect<T_,NC_> &x, LocalVect<T_,NR_> &y);
00220 
00225     void Mult(const LocalVect<T_,NC_> &x, LocalVect<T_,NR_> &y);
00226 
00229     void Symmetrize();
00230 
00239     int Factor();
00240 
00250     int Solve(LocalVect<T_,NR_> &b);
00251 
00261     int FactorAndSolve(LocalVect<T_,NR_> &b);
00262 
00265     void Invert(LocalMatrix<T_,NR_,NC_> &A);
00266 
00273     T_ getInnerProduct(const LocalVect<T_,NC_> &x, const LocalVect<T_,NR_> &y);
00274 
00276     T_ *getArray() { return _a; }
00277 
00278  private :
00279 
00280     mutable ErrorInLocalMatrix<T_,NR_,NC_> _e;
00281     size_t _length;
00282     T_ _a[NR_*NC_];
00283 };
00284 
00285 
00287 //                         I M P L E M E N T A T I O N                       //
00289 
00290 template<class T_,size_t NR_,size_t NC_>
00291 LocalMatrix<T_,NR_,NC_>::LocalMatrix()
00292 {
00293    _length = NR_ * NC_;
00294    for (size_t k=0; k<_length; k++)
00295       _a[k] = 0;
00296 }
00297 
00298 
00299 template<class T_,size_t NR_,size_t NC_>
00300 LocalMatrix<T_,NR_,NC_>::LocalMatrix(const LocalMatrix<T_,NR_,NC_> &m)
00301 {
00302    _length = m._length;
00303    for (size_t k=0; k<_length; k++)
00304       _a[k] = m._a[k];
00305 }
00306 
00307 
00308 template<class T_,size_t NR_,size_t NC_>
00309 LocalMatrix<T_,NR_,NC_>::LocalMatrix(Element *el, const SpMatrix<T_> &a)
00310 {
00311    Localize(el,a);
00312 #ifdef _OFELI_DEBUG
00313    std::clog << "An instance of class LocalMatrix is constructed.\n";
00314    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00315 #endif
00316 }
00317 
00318 
00319 template<class T_,size_t NR_,size_t NC_>
00320 LocalMatrix<T_,NR_,NC_>::LocalMatrix(Element *el, const SkMatrix<T_> &a)
00321 {
00322    Localize(el,a);
00323 #ifdef _OFELI_DEBUG
00324    std::clog << "An instance of class LocalMatrix is constructed.\n";
00325    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00326 #endif
00327 }
00328 
00329 
00330 template<class T_,size_t NR_,size_t NC_>
00331 LocalMatrix<T_,NR_,NC_>::LocalMatrix(Element *el, const SkSMatrix<T_> &a)
00332 {
00333    Localize(el,a);
00334 #ifdef _OFELI_DEBUG
00335    std::clog << "An instance of class NodeVect is constructed.\n";
00336    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00337 #endif
00338 }
00339 
00340 
00341 template<class T_,size_t NR_,size_t NC_>
00342 LocalMatrix<T_,NR_,NC_>::~LocalMatrix()
00343 {
00344 #ifdef _OFELI_DEBUG
00345    std::clog << "An instance of class LocalMatrix is destructed.\n";
00346    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00347 #endif
00348 }
00349 
00350 
00351 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00352 template<class T_,size_t NR_,size_t NC_>
00353 void LocalMatrix<T_,NR_,NC_>::set(int opt)
00354 {
00355    if (opt==IDENTITY) {
00356       if (NR_!=NC_)
00357          _e.Message(__FILE__,__LINE__,31);
00358       for (size_t i=1; i<=NR_*NC_; i++)
00359          _a[i] = 0;
00360       for (size_t j=1; j<=NR_; j++)
00361          _a[(j-1)*(NC_+1)] = 1;
00362    }
00363 }
00364 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00365 
00366 
00367 template<class T_,size_t NR_,size_t NC_>
00368 void LocalMatrix<T_,NR_,NC_>::Localize(Element *el, const SpMatrix<T_> &a)
00369 {
00370    size_t k = 0;
00371    for (size_t n=1; n<=el->getNbNodes(); n++) {
00372       Node *nd = el->getPtrNode(n);
00373       for (size_t i=1; i<=nd->getNbDOF(); i++)
00374          for (size_t j=1; j<=nd->getNbDOF(); j++)
00375             (*this)(i,j) = a(nd->getDOF(i),nd->getDOF(j));
00376    }
00377 }
00378 
00379 
00380 template<class T_,size_t NR_,size_t NC_>
00381 void LocalMatrix<T_,NR_,NC_>::Localize(Element *el, const SkMatrix<T_> &a)
00382 {
00383    size_t k = 0;
00384    for (size_t n=1; n<=el->getNbNodes(); n++) {
00385       Node *nd = el->getPtrNode(n);
00386       for (size_t i=1; i<=nd->getNbDOF(); i++)
00387          for (size_t j=1; j<=nd->getNbDOF(); j++)
00388             (*this)(i,j) = a(nd->getDOF(i),nd->getDOF(j));
00389    }
00390 }
00391 
00392 
00393 template<class T_,size_t NR_,size_t NC_>
00394 void LocalMatrix<T_,NR_,NC_>::Localize(Element *el, const SkSMatrix<T_> &a)
00395 {
00396    size_t k = 0;
00397    for (size_t n=1; n<=el->getNbNodes(); n++) {
00398       Node *nd = el->getPtrNode(n);
00399       for (size_t i=1; i<=nd->getNbDOF(); i++)
00400          for (size_t j=1; j<=nd->getNbDOF(); j++)
00401             (*this)(i,j) = a(nd->getDOF(i),nd->getDOF(j));
00402    }
00403 }
00404 
00405 
00406 template<class T_,size_t NR_,size_t NC_>
00407 LocalMatrix<T_,NR_,NC_> & LocalMatrix<T_,NR_,NC_>::operator=(const LocalMatrix<T_,NR_,NC_> &m)
00408 {
00409    _length = m._length;
00410    for (size_t k=0; k<_length; k++)
00411       _a[k] = m._a[k];
00412    return *this;
00413 }
00414 
00415 
00416 template<class T_,size_t NR_,size_t NC_>
00417 LocalMatrix<T_,NR_,NC_> & LocalMatrix<T_,NR_,NC_>::operator=(const T_ &x)
00418 {
00419    for (size_t k=0; k<_length; k++)
00420       _a[k] = x;
00421    return *this;
00422 }
00423 
00424 
00425 template<class T_,size_t NR_,size_t NC_>
00426 LocalMatrix<T_,NR_,NC_> & LocalMatrix<T_,NR_,NC_>::operator+=(const LocalMatrix<T_,NR_,NC_> &m)
00427 {
00428    for (size_t k=0; k<_length; k++)
00429       _a[k] += m._a[k];
00430    return *this;
00431 }
00432 
00433 
00434 template<class T_,size_t NR_,size_t NC_>
00435 LocalMatrix<T_,NR_,NC_> & LocalMatrix<T_,NR_,NC_>::operator-=(const LocalMatrix<T_,NR_,NC_> &m)
00436 {
00437    for (size_t k=0; k<_length; k++)
00438       _a[k] -= m._a[k];
00439    return *this;
00440 }
00441 
00442 
00443 template<class T_,size_t NR_,size_t NC_>
00444 LocalVect<T_,NR_> LocalMatrix<T_,NR_,NC_>::operator*(LocalVect<T_,NC_> &x)
00445 {
00446    LocalVect<T_,NR_> v;
00447    size_t k=0;
00448    for (size_t i=0; i<NR_; i++)
00449       for (size_t j=0; j<NC_; j++)
00450          v[i] += _a[k++]*x[j];
00451    return v;
00452 }
00453 
00454 
00455 template<class T_,size_t NR_,size_t NC_>
00456 LocalMatrix<T_,NR_,NC_> & LocalMatrix<T_,NR_,NC_>::operator+=(const T_ &x)
00457 {
00458    for (size_t k=0; k<_length; k++)
00459       _a[k] += x;
00460    return *this;
00461 }
00462 
00463 
00464 template<class T_,size_t NR_,size_t NC_>
00465 LocalMatrix<T_,NR_,NC_> & LocalMatrix<T_,NR_,NC_>::operator-=(const T_ &x)
00466 {
00467    for (size_t k=0; k<_length; k++)
00468       _a[k] -= x;
00469    return *this;
00470 }
00471 
00472 
00473 template<class T_,size_t NR_,size_t NC_>
00474 LocalMatrix<T_,NR_,NC_> & LocalMatrix<T_,NR_,NC_>::operator*=(const T_ &x)
00475 {
00476    for (size_t k=0; k<_length; k++)
00477       _a[k] *= x;
00478    return *this;
00479 }
00480 
00481 
00482 template<class T_,size_t NR_,size_t NC_>
00483 LocalMatrix<T_,NR_,NC_> & LocalMatrix<T_,NR_,NC_>::operator/=(const T_ &x)
00484 {
00485    for (size_t k=0; k<_length; k++)
00486       _a[k] /= x;
00487    return *this;
00488 }
00489 
00490 
00491 template<class T_,size_t NR_,size_t NC_>
00492 void LocalMatrix<T_,NR_,NC_>::MultAdd(const LocalVect<T_,NC_> &x, LocalVect<T_,NR_> &y)
00493 {
00494    size_t k=0;
00495    for (size_t i=0; i<NR_; i++)
00496       for (size_t j=0; j<NC_; j++)
00497          y[i] += _a[k++] * x[j];
00498 }
00499 
00500 
00501 template<class T_,size_t NR_,size_t NC_>
00502 void LocalMatrix<T_,NR_,NC_>::MultAddScal(const T_ &a, const LocalVect<T_,NC_> &x, LocalVect<T_,NR_> &y)
00503 {
00504    size_t k=0;
00505    for (size_t i=0; i<NR_; i++)
00506       for (size_t j=0; j<NC_; j++)
00507          y[i] += a * _a[k++] * x[j];
00508 }
00509 
00510 
00511 template<class T_,size_t NR_,size_t NC_>
00512 void LocalMatrix<T_,NR_,NC_>::Mult(const LocalVect<T_,NC_> &x, LocalVect<T_,NR_> &y)
00513 {
00514    y = 0;
00515    MultAdd(x,y);
00516 }
00517 
00518 
00519 template<class T_,size_t NR_,size_t NC_>
00520 void LocalMatrix<T_,NR_,NC_>::Symmetrize()
00521 {
00522    for (size_t i=0; i<NR_; i++)
00523       for (size_t j=0; j<i; j++)
00524          _a[i*NC_+j] = _a[j*NC_+i];
00525 }
00526 
00527 
00528 template<class T_,size_t NR_,size_t NC_>
00529 int LocalMatrix<T_,NR_,NC_>::Factor()
00530 {
00531    register size_t j, k;
00532    if (NR_!=NC_)
00533       _e.Message(__FILE__,__LINE__,11);
00534    for (size_t i=1; i<NR_; ++i) {
00535       for (j=1; j<=i; j++) {
00536          if (Abs(_a[NR_*(j-1)+j-1]) < OFELI_EPSMCH)
00537             _e.Message(__FILE__,__LINE__,12,j,0,_a[NR_*(j-1)+j-1]);
00538          _a[NR_*i+j-1] /= _a[NR_*(j-1)+j-1];
00539          for (k=0; k<j; k++)
00540             _a[NR_*i+j] -= _a[NR_*i+k]*_a[NR_*k+j];
00541       }
00542       for (j=i+1; j<NR_; ++j)
00543          for (k=0; k<i; ++k)
00544             _a[NR_*i+j] -= _a[NR_*i+k]*_a[NR_*k+j];
00545    }
00546    return 0;
00547 }
00548 
00549 
00550 template<class T_,size_t NR_,size_t NC_>
00551 int LocalMatrix<T_,NR_,NC_>::Solve(LocalVect<T_,NR_> &b)
00552 {
00553    register int i, j;
00554    if (NR_!=NC_)
00555       _e.Message(__FILE__,__LINE__,22);
00556    for (i=0; i<int(NR_); i++) {
00557       T_ s = 0;
00558       for (j=0; j<i; j++)
00559          s += _a[NR_*i+j] * b[j];
00560       b[i] -= s;
00561    }
00562    for (i=NR_-1; i>-1; i--) {
00563       if (Abs(_a[NR_*i+i]) < OFELI_EPSMCH)
00564          _e.Message(__FILE__,__LINE__,21,i+1,0,_a[NR_*i+i]);
00565       b[i] /= _a[NR_*i+i];
00566       for (j=0; j<i; j++)
00567          b[j] -= b[i] * _a[NR_*j+i];
00568    }
00569    return 0;
00570 }
00571 
00572 
00573 template<class T_,size_t NR_,size_t NC_>
00574 int LocalMatrix<T_,NR_,NC_>::FactorAndSolve(LocalVect<T_,NR_> &b)
00575 {
00576    Factor();
00577    Solve(b);
00578    return 0;
00579 }
00580 
00581 
00582 template<class T_,size_t NR_,size_t NC_>
00583 void LocalMatrix<T_,NR_,NC_>::Invert(LocalMatrix<T_,NR_,NC_> &A)
00584 {
00585    LocalVect<T_,NR_> b;
00586    Factor();
00587    for (size_t i=1; i<=NR_; i++) {
00588       b = 0;
00589       b(i) = 1;
00590       Solve(b);
00591       for (size_t j=1; j<=NR_; j++)
00592          A(j,i) = b(j);
00593    }
00594 }
00595 
00596 
00597 template<class T_,size_t NR_,size_t NC_>
00598 T_ LocalMatrix<T_,NR_,NC_>::getInnerProduct(const LocalVect<T_,NC_> &x, const LocalVect<T_,NR_> &y)
00599 {
00600    double s = 0;
00601    for (size_t i=1; i<=NC_; i++)
00602       for (size_t j=1; i<=NR_; i++)
00603          s += _a[(i-1)*NC_+j-1]*x(i)*y(j);
00604    return s;
00605 }
00606 
00607 
00609 //              A S S O C I A T E D    F U N C T I O N S                     //
00611 
00617 template<class T_,size_t NR_,size_t NC_>
00618 LocalMatrix<T_,NR_,NC_> operator*(T_ a, const LocalMatrix<T_,NR_,NC_> &x)
00619 {
00620    LocalMatrix<T_,NR_,NC_> z;
00621    for (size_t i=1; i<=NR_; ++i)
00622       for (size_t j=1; j<=NC_; ++j)
00623          z(i,j) = a*x(i,j);
00624    return z;
00625 }
00626 
00627 
00633 template<class T_,size_t NR_,size_t NC_>
00634 LocalMatrix<T_,NR_,NC_> operator/(T_ a, const LocalMatrix<T_,NR_,NC_> &x)
00635 {
00636    LocalMatrix<T_,NR_,NC_> z;
00637    for (size_t i=1; i<=NR_; ++i)
00638       for (size_t j=1; j<=NC_; ++j)
00639          z(i,j) = x(i,j)/a;
00640    return z;
00641 }
00642 
00643 
00649 template<class T_,size_t NR_,size_t NC_>
00650 LocalMatrix<T_,NR_,NC_> operator+(const LocalMatrix<T_,NR_,NC_> &x, const LocalMatrix<T_,NR_,NC_> &y)
00651 {
00652    LocalMatrix<T_,NR_,NC_> z;
00653    for (size_t i=1; i<=NR_; ++i)
00654       for (size_t j=1; j<=NC_; ++j)
00655          z(i,j) = x(i,j) + y(i,j);
00656    return z;
00657 }
00658 
00659 
00665 template<class T_,size_t NR_,size_t NC_>
00666 LocalMatrix<T_,NR_,NC_> operator-(const LocalMatrix<T_,NR_,NC_> &x, const LocalMatrix<T_,NR_,NC_> &y)
00667 {
00668    LocalMatrix<T_,NR_,NC_> z;
00669    for (size_t i=0; i<NR_; ++i)
00670       for (size_t j=1; j<=NC_; ++j)
00671          z(i,j) = x(i,j) - y(i,j);
00672    return z;
00673 }
00674 
00675 
00679 template<class T_, size_t NR_, size_t NC_>
00680 ostream& operator<<(ostream &s, const LocalMatrix<T_,NR_,NC_> &a)
00681 {
00682    s.width(6);
00683    s.setf(ios::right|ios::scientific);
00684    s << endl;
00685    for (size_t i=1; i<=NR_; i++) {
00686       s << "\nRow  " << setw(6) << i << endl;
00687       for (size_t j=1; j<=NC_; j++)
00688          s << "  " << setprecision(8) << std::setfill(' ') << setw(18) << a(i,j);
00689       s << endl;
00690    }
00691    return s;
00692 }
00693 
00694 } /* namespace OFELI */
00695 
00696 #endif