DSMatrix.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 class DSMatrix for symmetric matrices
00030 
00031   ==============================================================================*/
00032 
00033 
00034 #ifndef __DSMATRIX_H
00035 #define __DSMATRIX_H
00036 
00037 #include "Matrix.h"
00038 
00039 
00040 namespace OFELI {
00041 
00042 template<class T_> class DSMatrix;
00043 
00048 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00049 template<class T_>
00050 struct ErrorInDSMatrix {
00051    void Message(const char *file, size_t line, int code, int p1=0, int p2=0, T_ p3=T_(0)) {
00052      cerr << "\n*** Fatal Error in OFELI ***\n";
00053      cerr << "File : " << file << ", line : " << line << "\nIn DMatrix<>::";
00054      switch (code) {
00055 
00056         case 31: 
00057            cerr << "Factor() : The " << p1 << "-th pivot = " << p3 << " is too small\n";
00058            break;
00059 
00060         case 72: 
00061            cerr << "Solve(b) : Warning : Matrix has not been factorized\n";
00062            return;
00063      }
00064      cerr << "Program stops" << endl;
00065      exit(1);
00066    }
00067 };
00068 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00069 
00070 
00083 template<class T_>
00084 class DSMatrix : public Matrix<T_>
00085 {
00086 
00087   using Matrix<T_>::_rmin;
00088   using Matrix<T_>::_cmin;
00089   using Matrix<T_>::_rmax;
00090   using Matrix<T_>::_cmax;
00091   using Matrix<T_>::_nb_rows;
00092   using Matrix<T_>::_nb_cols;
00093   using Matrix<T_>::_size;
00094   using Matrix<T_>::_length;
00095   using Matrix<T_>::_a;
00096   using Matrix<T_>::_diag;
00097   using Matrix<T_>::_zero;
00098   using Matrix<T_>::_ch;
00099   using Matrix<T_>::_fact;
00100 
00101  public :
00102 
00105     DSMatrix();
00106 
00111     DSMatrix(size_t dim);
00112 
00115     DSMatrix(const DSMatrix<T_> &m);
00116 
00118     ~DSMatrix();
00119 
00121     void setDiag();
00122 
00125     void Resize(size_t dim);
00126 
00132     void Set(size_t i, size_t j, const T_ &x);
00133 
00139     void Add(size_t i, size_t j, const T_ &x);
00140 
00145     T_ operator()(size_t i, size_t j) const;
00146 
00151     T_ & operator()(size_t i, size_t j);
00152 
00155     DSMatrix<T_> & operator=(const DSMatrix<T_> &m);
00156 
00159     DSMatrix<T_> & operator=(const T_ &x);
00160 
00168     int Factor();
00169 
00171     void MultAdd(const Vect<T_> &x, Vect<T_> &y) const;
00172 
00178     void MultAdd(T_ a, const Vect<T_> &x, Vect<T_> &y) const;
00179 
00181     void Mult(const Vect<T_> &x, Vect<T_> &y) const;
00182 
00191     int Solve(Vect<T_> &b);
00192 
00203     int Solve(const Vect<T_> &b, Vect<T_> &x);
00204 
00205 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00207     size_t getColInd(size_t i) const;
00208 
00210     size_t getRowPtr(size_t i) const;
00211 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00212 
00217     T_ *getArray() const { return _a; }
00218 
00220     T_ getEntry(size_t i, size_t j) const;
00221 
00228     void setPrintView(size_t rmin, size_t rmax, size_t cmin, size_t cmax);
00229 
00230  protected:
00231 
00232 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00233    mutable ErrorInDSMatrix<T_> _e;
00234    void set_(size_t i, size_t j, T_ x) { _a[(i-1)*(i-2)/2+j-1] = x; }
00235 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00236 
00237 };
00238 
00240 //                         I M P L E M E N T A T I O N                       //
00242 
00243 
00244 template<class T_>
00245 DSMatrix<T_>::DSMatrix()
00246 {
00247    _nb_rows = _length = _size = 0;
00248    _a = NULL;
00249    _fact = false;
00250 #ifdef _OFELI_DEBUG
00251    std::clog << "An instance of class DSMatrix is constructed.\n";
00252    std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
00253 #endif
00254 }
00255 
00256 
00257 template<class T_>
00258 DSMatrix<T_>::DSMatrix(size_t dim)
00259 {
00260    Resize(dim);
00261    _rmin = 1;
00262    _rmax = dim;
00263    _cmin = 1;
00264    _cmax = dim;
00265 #ifdef _OFELI_DEBUG
00266    std::clog << "An instance of class DSMatrix is constructed.\n";
00267    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00268 #endif
00269 }
00270 
00271 
00272 template<class T_>
00273 DSMatrix<T_>::DSMatrix(const DSMatrix<T_> &m)
00274 {
00275    _size = _nb_rows = _nb_cols = m._size;
00276    _length = m._length;
00277    _a.resize(_length);
00278    _a = m._a;
00279    _rmin = m._rmin;
00280    _rmax = m._rmax;
00281    _cmin = m._cmin;
00282    _cmax = m._cmax;
00283 #ifdef _OFELI_DEBUG
00284    std::clog << "An instance of class DSMatrix is constructed.\n";
00285    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00286 #endif
00287 }
00288 
00289 
00290 template<class T_>
00291 DSMatrix<T_>::~DSMatrix()
00292 {
00293 #ifdef _OFELI_DEBUG
00294    std::clog << "An instance of class DSMatrix is destructed.\n";
00295    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00296 #endif
00297 }
00298 
00299 
00300 template<class T_>
00301 void DSMatrix<T_>::setDiag()
00302 {
00303    for (size_t i=0; i<_size; i++)
00304       _diag[i] = _a[_nb_cols*i+i];
00305 }
00306 
00307 
00308 template<class T_>
00309 void DSMatrix<T_>::Resize(size_t dim)
00310 {
00311    _size = _nb_rows = _nb_cols = dim;
00312    _length = _size*(_size+1)/2;
00313    _a.resize(_length);
00314    Clear(_a);
00315 }
00316 
00317 
00318 template<class T_>
00319 void DSMatrix<T_>::Set(size_t i, size_t j, const T_ &x)
00320 {
00321 #ifdef _BOUNDS_DEBUG
00322    assert(i>0);
00323    assert(i<=_nb_rows);
00324    assert(j>0);
00325    assert(j<=_nb_rows);
00326 #endif
00327    if (i>=j)
00328       _a[i*(i-1)/2+j-1] = x;
00329 }
00330 
00331 
00332 template<class T_>
00333 void DSMatrix<T_>::Add(size_t i, size_t j, const T_ &x)
00334 {
00335 #ifdef _BOUNDS_DEBUG
00336    assert(i>0);
00337    assert(i<=_nb_rows);
00338    assert(j>0);
00339    assert(j<=_nb_rows);
00340 #endif
00341    if (i>=j)
00342       _a[i*(i-1)/2+j-1] += x;
00343 }
00344 
00345 
00346 template<class T_>
00347 T_ DSMatrix<T_>::operator()(size_t i, size_t j) const
00348 {
00349 #ifdef _BOUNDS_DEBUG
00350    assert(i>0);
00351    assert(i<=_nb_rows);
00352    assert(j>0);
00353    assert(j<=_nb_rows);
00354 #endif
00355    if (i<j)
00356       return _a[j*(j-1)/2+i-1];
00357    else
00358       return _a[i*(i-1)/2+j-1];
00359 }
00360 
00361 
00362 template<class T_>
00363 T_ & DSMatrix<T_>::operator()(size_t i, size_t j)
00364 {
00365 #ifdef _BOUNDS_DEBUG
00366    assert(i>0);
00367    assert(i<=_nb_rows);
00368    assert(j>0);
00369    assert(j<=_nb_cols);
00370 #endif
00371    if (i>=j)
00372       return _a[i*(i-1)/2+j-1];
00373    else
00374       return _zero;
00375 }
00376 
00377 
00378 template<class T_>
00379 DSMatrix<T_> & DSMatrix<T_>::operator=(const DSMatrix<T_> &m)
00380 {
00381    _a = m._a;
00382    return *this;
00383 }
00384 
00385 
00386 template<class T_>
00387 DSMatrix<T_> & DSMatrix<T_>::operator=(const T_ &x)
00388 {
00389    for (size_t i=0; i<_length; i++)
00390       _a[i] = x;
00391    return *this;
00392 }
00393 
00394 
00395 template<class T_>
00396 int DSMatrix<T_>::Factor()
00397 {
00398    register size_t i, j, k;
00399    int err = 0;
00400    T_ s, pivot;
00401    for (i=0; i<_size; i++) {
00402       for (j=0; j<i; j++)
00403          for (k=0; k<j; k++)
00404             _a[(i+1)*i/2+j] -= _a[(i+1)*i/2+k]*_a[(j+1)*j/2+k];
00405       pivot = _a[(i+1)*i/2+i];
00406       for (j=0; j<i; j++) {
00407          s = _a[(i+1)*i/2+j]*_a[(j+1)*j/2+j];
00408          pivot -= s*_a[(i+1)*i/2+j];
00409          _a[(i+1)*i/2+j] = s;
00410       }
00411       if (Abs(pivot) < OFELI_TOLERANCE)
00412          _e.Message(__FILE__,__LINE__,31,i+1,0,pivot);
00413       if (pivot < 0)
00414          err = -1;
00415       _a[(i+1)*i/2+i] = 1./pivot;
00416    }
00417    _fact = true;
00418    return err;
00419 }
00420 
00421 
00422 template<class T_>
00423 void DSMatrix<T_>::MultAdd(const Vect<T_> &x, Vect<T_> &y) const
00424 {
00425    for (size_t i=1; i<=_nb_rows; i++) {
00426       for (size_t j=1; j<=i; j++)
00427          y(i) += _a[(i-1)*i/2+j-1] * x(j);
00428       for (size_t k=i+1; k<=_nb_rows; k++)
00429          y(i) += _a[(k-1)*k/2+i-1] * x(k);
00430    }
00431 }
00432 
00433 
00434 template<class T_>
00435 void DSMatrix<T_>::MultAdd(T_ a, const Vect<T_> &x, Vect<T_> &y) const
00436 {
00437    for (size_t i=1; i<=_nb_rows; i++) {
00438       for (size_t j=1; j<=i; j++)
00439          y(i) += a * _a[(i-1)*i/2+j-1] * x(j);
00440       for (size_t k=i+1; k<=_nb_rows; k++)
00441          y(i) += a * _a[(k-1)*k/2+i-1] * x(k);
00442    }
00443 }
00444 
00445 
00446 template<class T_>
00447 void DSMatrix<T_>::Mult(const Vect<T_> &x, Vect<T_> &y) const
00448 {
00449    y = 0;
00450    MultAdd(x,y);
00451 }
00452 
00453 
00454 template<class T_>
00455 int DSMatrix<T_>::Solve(Vect<T_> &b)
00456 {
00457    if (!_fact)
00458       _e.Message(__FILE__,__LINE__,72);
00459    int i, j;
00460    for (i=0; i<_size; i++) {
00461       T_ s = 0;
00462       for (j=0; j<i; j++)
00463          s += _a[(i+1)*i/2+j] * b[j];
00464       b[i] -= s;
00465    }
00466    for (i=0; i<_size; i++)
00467       b[i] *= _a[(i+1)*i/2+i];
00468    for (i=_size-1; i>-1; i--)
00469       for (j=0; j<i; j++)
00470          b[j] -= b[i] * _a[(i+1)*i/2+j];
00471    return 0;
00472 }
00473 
00474 
00475 template<class T_>
00476 int DSMatrix<T_>::Solve(const Vect<T_> &b, Vect<T_> &x)
00477 {
00478    x = b;
00479    return Solve(x);
00480 }
00481 
00482 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00483 template<class T_>
00484 size_t DSMatrix<T_>::getColInd(size_t i) const
00485 {
00486 #ifdef _BOUNDS_DEBUG
00487    assert(i<=_length && i>0);
00488 #endif
00489    return _ch[i-1];
00490 }
00491 
00492 
00493 template<class T_>
00494 size_t DSMatrix<T_>::getRowPtr(size_t i) const
00495 {
00496 #ifdef _BOUNDS_DEBUG
00497    assert(i<=_nb_rows+1 && i>0);
00498 #endif
00499    return _ch[i-1];
00500 }
00501 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00502 
00503 
00504 template<class T_>
00505 T_ DSMatrix<T_>::getEntry(size_t i, size_t j) const
00506 {
00507    if (i>=j)
00508       return _a[i*(i-1)/2+j-1];
00509    else
00510       return _a[j*(j-1)/2+i-1];
00511 }
00512 
00513 
00514 template<class T_>
00515 void DSMatrix<T_>::setPrintView(size_t rmin, size_t rmax, size_t cmin, size_t cmax)
00516 {
00517    _rmin = rmin;
00518    _rmax = rmax;
00519    if (_rmin==0)
00520       _rmin = 1;
00521    if (_rmin>_size)
00522       _rmin = _size;
00523    _cmin = cmin;
00524    _cmax = cmax;
00525    if (_cmin==0)
00526       _cmin = 1;
00527    if (_cmin>_size)
00528       _cmin = _size;
00529 }
00530 
00532 //                 A S S O C I A T E D   F U N C T I O N S                   //
00534 
00535 
00539 template<class T_>
00540 ostream& operator<<(ostream& s, const DSMatrix<T_> & a)
00541 {
00542    size_t rmin, rmax, cmin, cmax;
00543    a.getPrintView(rmin,rmax,cmin,cmax);
00544    s.setf(ios::right|ios::scientific);
00545    s << endl;
00546    for (size_t i=rmin-1; i<rmax; i++) {
00547       s << "\nRow  " << setw(6) << i+1 << endl;
00548       for (size_t j=cmin-1; j<=cmax; j++)
00549          s << "  " << setprecision(8) << std::setfill(' ') << setw(18) << a.Array()[i*(i+1)/2+j];
00550       s << endl;
00551    }
00552    return s;
00553 }
00554 
00555 } /* namespace OFELI */
00556 
00557 #endif