TrMatrix.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                    Template Class Matrix for tridiagonal matrices
00030 
00031   ==============================================================================*/
00032 
00033 
00034 #ifndef _TR_MATRIX_H
00035 #define _TR_MATRIX_H
00036 
00037 
00038 #include "Matrix.h"
00039 
00040 
00041 namespace OFELI {
00042 
00047 template<class T_> class TrMatrix;
00048 
00049 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00050 template<class T_>
00051 struct ErrorInTrMatrix {
00052   void Message(const char *file, size_t line, int code, int p1=0, int p2=0, double p3=0) {
00053      cerr << "\n*** Fatal Error in OFELI ***\n";
00054      cerr << "File : " << file << ", line : " << line << "\nIn TrMatrix<>::";
00055      switch (code) {
00056 
00057         case 31: 
00058            cerr << "Factor(): matrix  cannot be fatorized separately.\n";
00059            cerr << "Call Solve(b) directly.\n";
00060            return;
00061 
00062         case 41: 
00063            cerr << "Solve(b): The " << p1 << "-th pivot = " << p3 << " is too small\n";
00064            break;
00065 
00066         case 51:
00067            cerr << "Operator (): Index pair (" << p1 << "," << p2 << ") is not compatible with tridiagonal structure\n";
00068            break;
00069      }
00070      cerr << "Program stops" << endl;
00071      exit(1);
00072   }
00073 };
00074 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00075 
00076 
00090 template<class T_>
00091 class TrMatrix : public Matrix<T_>
00092 {
00093 
00094  public :
00095 
00097     TrMatrix();
00098 
00101     TrMatrix(size_t size);
00102 
00104     TrMatrix(const TrMatrix &m);
00105 
00107     ~TrMatrix();
00108 
00109 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00110     virtual void setMesh(const class Mesh &mesh, size_t dof=0, size_t type=0);
00111 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00112 
00113 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00114     virtual void setMesh(size_t dof, const class Mesh &mesh, int code=0);
00115 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00116 
00117 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00118     virtual void setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh);
00119 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00120 
00121 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00122     void setDiag() { }
00123 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00124 
00127     void Resize(size_t size);
00128 
00130     void MultAdd(const Vect<T_> &x, Vect<T_> &y) const;
00131 
00133     void MultAdd(T_ a, const Vect<T_> &x, Vect<T_> &y) const;
00134 
00136     void Mult(const Vect<T_> &x, Vect<T_> &y) const;
00137 
00139     void TMult(const Vect<T_> &x, Vect<T_> &y) const;
00140 
00142     void Set(size_t i, size_t j, const T_ &x);
00143 
00145     void Add(size_t i, size_t j, const T_ &x);
00146 
00151     const T_ operator()(size_t i) const { return _a[i-1]; }
00152 
00156     T_ operator()(size_t i, size_t j) const;
00157 
00162     T_ & operator()(size_t i, size_t j);
00163 
00166     TrMatrix<T_> & operator=(const TrMatrix<T_> &m);
00167 
00170     TrMatrix<T_> & operator=(const T_ &x);
00171 
00174     TrMatrix<T_> & operator*=(const T_ &x);
00175 
00176 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00177     int Factor();
00178 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00179 
00190     int Solve(Vect<T_> &b);
00191 
00192 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00194     size_t getColInd(size_t i) const;
00195 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00196 
00197 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00199     size_t getRowPtr(size_t i) const;
00200 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00201 
00204     T_ *getArray() const { return _a; }
00205 
00207     T_ getEntry(size_t i, size_t j) const;
00208 
00209  private :
00210 
00211    size_t _imin, _imax;
00212    mutable ErrorInTrMatrix<T_> _e;
00213 
00214    using Matrix<T_>::_rmin;
00215    using Matrix<T_>::_cmin;
00216    using Matrix<T_>::_rmax;
00217    using Matrix<T_>::_cmax;
00218    using Matrix<T_>::_size;
00219    using Matrix<T_>::_nb_cols;
00220    using Matrix<T_>::_nb_rows;
00221    using Matrix<T_>::_length;
00222    using Matrix<T_>::_a;
00223    using Matrix<T_>::_diag;
00224    using Matrix<T_>::_fact;
00225    using Matrix<T_>::_zero;
00226    using Matrix<T_>::_temp;
00227    using Matrix<T_>::_ch;
00228 };
00229 
00231 //                         I M P L E M E N T A T I O N                       //
00233 
00234 template<class T_>
00235 TrMatrix<T_>::TrMatrix()
00236 {
00237    _size = 0;
00238 #ifdef _OFELI_DEBUG
00239    std::clog << "An instance of class TrMatrix is constructed.\n";
00240    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00241 #endif
00242 }
00243 
00244 template<class T_>
00245 TrMatrix<T_>::TrMatrix(size_t size)
00246 {
00247    _size = _nb_rows = _nb_cols = size;
00248    _length = 3*_size;
00249    _a.resize(_length);
00250    Clear(_a);
00251    _ch.resize(_size);
00252    _ch[0] = 0;
00253    for (size_t i=1; i<_size; i++)
00254       _ch[i] = i+1;
00255 }
00256 
00257 template<class T_>
00258 TrMatrix<T_>::TrMatrix(const TrMatrix &m)
00259 {
00260    _size = m._size;
00261    _length = m._length;
00262    _a.resize(_length);
00263    _a = m._a;
00264 }
00265 
00266 template<class T_>
00267 TrMatrix<T_>::~TrMatrix()
00268 {
00269 #ifdef _OFELI_DEBUG
00270    std::clog << "An instance of class Trmatrix is destructed.\n";
00271    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00272 #endif
00273 }
00274 
00275 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00276 template<class T_>
00277 void TrMatrix<T_>::setMesh(const class Mesh &mesh, size_t dof, size_t type)
00278 {
00279    if (mesh.getDim()==0) { }
00280    type = 0;
00281    dof = 0;
00282    _length = 3*_size;
00283    _diag.resize(_size);
00284    _a.resize(_length);
00285    Clear(_a);
00286    _zero = T_(0);
00287    _fact = false;
00288 }
00289 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00290 
00291 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00292 template<class T_>
00293 void TrMatrix<T_>::setMesh(size_t dof, const class Mesh &mesh, int code)
00294 {
00295 // This is just to avoid warning on unused variable
00296    dof  = 0;
00297    code = 0;
00298    if (mesh.getDim()==0) { }
00299 }
00300 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00301 
00302 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00303 template<class T_>
00304 void TrMatrix<T_>::setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh)
00305 {
00306 // This is just to avoid warning on unused variable
00307    dof = 0;
00308    nb_eq = 0;
00309    if (mesh.getDim()==0) { }
00310 }
00311 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00312 
00313 template<class T_>
00314 void TrMatrix<T_>::Resize(size_t size)
00315 {
00316    _nb_rows = _nb_cols = _size = size;
00317    _length = 3 * _nb_rows;
00318    _a.resize(_length);
00319    Clear(_a);
00320 }
00321 
00322 
00323 template<class T_>
00324 void TrMatrix<T_>::MultAdd(const Vect<T_> &x, Vect<T_> &y) const
00325 {
00326    y(1) += _a[1]*x(1) + _a[2]*x(2);
00327    for (size_t i=2; i<_size; ++i)
00328       y(i) += _a[3*i-3]*x(i-1) + _a[3*i-2]*x(i) + _a[3*i-1]*x(i+1);
00329    y(_size) += _a[3*_size-3]*x(_size-1) + _a[3*_size-2]*x(_size);
00330 }
00331 
00332 
00333 template<class T_>
00334 void TrMatrix<T_>::MultAdd(T_ a, const Vect<T_> &x, Vect<T_> &y) const
00335 {
00336    y(1) += _a[1]*x(1) + _a[2]*x(2);
00337    for (size_t i=2; i<_size; ++i)
00338       y(i) += a * _a[3*i-3]*x(i-1) + _a[3*i-2]*x(i) + _a[3*i-1]*x(i+1);
00339    y(_size) += a * _a[3*_size-3]*x(_size-1) + _a[3*_size-2]*x(_size);
00340 }
00341 
00342 
00343 template<class T_>
00344 void TrMatrix<T_>::Mult(const Vect<T_> &x, Vect<T_> &y) const
00345 {
00346    y = 0;
00347    MultAdd(x,y);
00348 }
00349 
00350 
00351 template<class T_>
00352 void TrMatrix<T_>::TMult(const Vect<T_> &x, Vect<T_> &y) const
00353 {
00354    y(1) = _a[1]*x(1) + _a[3]*x(2);
00355    for (size_t i=2; i<_size; ++i)
00356       y(i) = _a[3*i-4]*x(i-1) + _a[3*i-2]*x(i) + _a[3*i]*x(i+1);
00357    y(_size) = _a[3*_size-4]*x(_size-1) + _a[3*_size-2]*x(_size);
00358 }
00359 
00360 
00361 template<class T_>
00362 void TrMatrix<T_>::Set(size_t i, size_t j, const T_ &x)
00363 {
00364    if (i==j)
00365       _a[3*i-2] = x;
00366    else if (i==j+1)
00367       _a[3*i-3] = x;
00368    else if (i==j-1)
00369       _a[3*i-1] = x;
00370 }
00371 
00372 
00373 template<class T_>
00374 void TrMatrix<T_>::Add(size_t i, size_t j, const T_ &x)
00375 {
00376    if (i==j)
00377       _a[3*i-2] += x;
00378    else if (i==j+1)
00379       _a[3*i-3] += x;
00380    else if (i==j-1)
00381       _a[3*i-1] += x;
00382 }
00383 
00384 
00385 template<class T_>
00386 T_ TrMatrix<T_>::operator()(size_t i, size_t j) const
00387 {
00388 #ifdef _BOUNDS_DEBUG
00389    assert(i==j || i==j+1 || i==j-1);
00390 #endif
00391    if (i==j)
00392       return _a[3*i-2];
00393    else if (i==j+1)
00394       return _a[3*i-3];
00395    else if (i==j-1)
00396       return _a[3*i-1];
00397    else
00398       return _zero;
00399 }
00400 
00401 
00402 template<class T_>
00403 T_ & TrMatrix<T_>::operator()(size_t i, size_t j)
00404 {
00405 #ifdef _BOUNDS_DEBUG
00406    assert(i==j || i==j+1 || i==j-1);
00407 #endif
00408    try {
00409       if (i==j)
00410          return _a[3*i-2];
00411       else if (i==j+1)
00412          return _a[3*i-3];
00413       else if (i==j-1)
00414          return _a[3*i-1];
00415       else
00416          throw(51);
00417    }
00418    catch(int n) { _e.Message(__FILE__,__LINE__,n,int(i),int(j)); }
00419    return _zero;
00420 }
00421 
00422 
00423 template<class T_>
00424 TrMatrix<T_> & TrMatrix<T_>::operator=(const TrMatrix<T_> &m)
00425 {
00426    _a = m._a;
00427    return *this;
00428 }
00429 
00430 
00431 template<class T_>
00432 TrMatrix<T_> & TrMatrix<T_>::operator=(const T_ &x)
00433 {
00434    for (size_t i=0; i<_length; ++i)
00435       _a[i] = x;
00436    return *this;
00437 }
00438 
00439 
00440 template<class T_>
00441 TrMatrix<T_> & TrMatrix<T_>::operator*=(const T_ &x)
00442 {
00443    for (size_t i=0; i<_length; i++)
00444       _a[i] *= x;
00445    return *this;
00446 }
00447 
00448 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00449 template<class T_>
00450 int TrMatrix<T_>::Factor()
00451 {
00452    _e.Message(__FILE__,__LINE__,31);
00453    return 1;
00454 }
00455 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00456 
00457 
00458 template<class T_>
00459 int TrMatrix<T_>::Solve(Vect<T_> &b)
00460 {
00461    for (size_t i=1; i<_size; ++i) {
00462       if (Abs(_a[3*i-2]) < OFELI_TOLERANCE)
00463          _e.Message(__FILE__,__LINE__,41,i+1,0,_a[3*i-2]);
00464       T_ p = _a[3*i]/_a [3*i-2];
00465       _a[3*i+1] -= p*_a[3*i-1];
00466       b(i+1) -= p*b(i);
00467    }
00468    if (Abs(_a[3*_size-2]) < OFELI_TOLERANCE)
00469       _e.Message(__FILE__,__LINE__,41,_size,0,_a[3*_size-2]);
00470    b(_size) /= _a[3*_size-2];
00471    for (int j=_size-1; j>0; j--)
00472       b(j) = (b(j)-_a[3*j-1]*b(j+1))/_a[3*j-2];
00473    return 0;
00474 }
00475 
00476 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00477 template<class T_>
00478 size_t TrMatrix<T_>::getColInd(size_t i) const
00479 {
00480 #ifdef _BOUNDS_DEBUG
00481    assert(i<=_length && i>0);
00482 #endif
00483    return _ch[i-1];
00484 }
00485 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00486 
00487 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00488 template<class T_>
00489 size_t TrMatrix<T_>::getRowPtr(size_t i) const
00490 {
00491 #ifdef _BOUNDS_DEBUG
00492    assert(i<=_nb_rows+1 && i>0);
00493 #endif
00494    return _ch[i-1];
00495 }
00496 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00497 
00498 
00499 template<class T_>
00500 T_ TrMatrix<T_>::getEntry(size_t i, size_t j) const
00501 {
00502    if (i==j)
00503       return _a[3*i-2];
00504    else if (i==j+1)
00505       return _a[3*i-3];
00506    else if (i==j-1)
00507       return _a[3*i-1];
00508    else
00509       return _zero;
00510 }
00511 
00513 //                           ASSOCIATED  FUNCTIONS                           //
00515 
00516 
00522 template<class T_>
00523 TrMatrix<T_> operator*(T_ a, const TrMatrix<T_> &A)
00524 {
00525    TrMatrix<T_> v(A);
00526    for (size_t i=0; i<A.getLength(); ++i)
00527       v[i] *= a;
00528    return v;
00529 }
00530 
00535 template<class T_>
00536 ostream& operator<<(ostream& s, const TrMatrix<T_> & a)
00537 {
00538    size_t rmin, rmax, cmin, cmax;
00539    a.getPrintView(rmin,rmax,cmin,cmax);
00540    s.setf(ios::scientific);
00541    s << endl;
00542    for (size_t i=rmin; i<=rmax; i++) {
00543       s << "\nRow  " << setw(6) << i << endl;
00544       for (size_t j=cmin; j<=cmax; j++) {
00545          if (j<i-1 || j>i+1)
00546             s << "  " << setprecision(8) << std::setfill(' ') << setw(18) << 0.;
00547          else if (j==i-1)
00548             s << "  " << setprecision(8) << std::setfill(' ') << setw(18) << a(i,j);
00549          else if (j==i)
00550             s << "  " << setprecision(8) << std::setfill(' ') << setw(18) << a(i,j);
00551          else if (j==i+1)
00552             s << "  " << setprecision(8) << std::setfill(' ') << setw(18) << a(i,j);
00553       }
00554       s << endl;
00555    }
00556    return s;
00557 }
00558 
00559 } /* namespace OFELI */
00560 
00561 #endif