SkMatrix.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 'SkMatrix' for general skyline matrix
00030 
00031   ==============================================================================*/
00032 
00033 
00034 #ifndef __SKMATRIX_H
00035 #define __SKMATRIX_H
00036 
00037 
00038 #include "Matrix.h"
00039 
00040 namespace OFELI {
00041 
00046 template<class T_> class SkMatrix;
00047 
00048 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00049 template<class T_>
00050 struct ErrorInSkMatrix
00051 {
00052    void Message(const char *file, size_t line, int code, int p1=0, int p2=0, T_ p3=0)
00053    {
00054      cerr << "\n*** Fatal Error in OFELI ***\n";
00055      cerr << "File : " << file << ", line : " << line << "\nIn SkMatrix<>::";
00056      switch (code) {
00057 
00058         case 51:
00059            cerr << "Set(i,j,x) : Illegal pair of indices : " << p1 << "," << p2 << endl;
00060            break;
00061 
00062         case 61:
00063            cerr << "Factor() : The " << p1 << "th pivot = " << p3 << " is too small\n";
00064            break;
00065 
00066         case 71:
00067            cerr << "Solve(b) : The " << p1 << "th diagonal entry = " << p3 << " is too small\n";
00068            break;
00069 
00070         case 72:
00071            cerr << "Solve(b) : Warning : Matrix has not been factorized\n";
00072            return;
00073 
00074         case 81:
00075            cerr << "Operator () : Index pair (" << p1 << "," << p2 << ") is not compatible with skyline structure\n";
00076 break;
00077      }
00078      cerr << "Program stops" << endl;
00079      exit(1);
00080    }
00081 };
00082 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00083 
00084 
00098 template<class T_>
00099 class SkMatrix : public Matrix<T_>
00100 {
00101    using Matrix<T_>::_rmin;
00102    using Matrix<T_>::_cmin;
00103    using Matrix<T_>::_rmax;
00104    using Matrix<T_>::_cmax;
00105    using Matrix<T_>::_size;
00106    using Matrix<T_>::_length;
00107    using Matrix<T_>::_zero;
00108    using Matrix<T_>::_temp;
00109    using Matrix<T_>::_fact;
00110    using Matrix<T_>::_a;
00111    using Matrix<T_>::_aU;
00112    using Matrix<T_>::_diag;
00113    using Matrix<T_>::_ch;
00114    using Matrix<T_>::_dof_type;
00115 
00116  public:
00117 
00119     SkMatrix();
00120 
00129     SkMatrix(const class Mesh &mesh, size_t dof=0);
00130 
00136     SkMatrix(size_t size);
00137 
00142     SkMatrix(const Vect<size_t> &ColHt);
00143 
00145     SkMatrix(const SkMatrix<T_> &m);
00146 
00148     ~SkMatrix();
00149 
00160     virtual void setMesh(const class Mesh &mesh, size_t dof=0, size_t type=0);
00161 
00162 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00163     void setMesh(size_t dof, const class Mesh &mesh, int code=0);
00164 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00165 
00166 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00167     void setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh);
00168 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00169 
00174     void setSkyline(const class Mesh &mesh);
00175 
00177     void setDiag();
00178 
00181     void setDOF(size_t i) { _dof = i; }
00182 
00188     void Set(size_t i, size_t j, const T_ &x);
00189 
00190 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00191     void SSet(size_t i, size_t j, const T_ &x);
00192 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00193 
00198     void MultAdd(const Vect<T_> &x, Vect<T_> &y) const;
00199 
00204     void TMultAdd(const Vect<T_> &x, Vect<T_> &y) const;
00205 
00211     void MultAdd(T_ a, const Vect<T_> &v, Vect<T_> &w) const;
00212 
00217     void Mult(const Vect<T_> &x, Vect<T_> &y) const;
00218 
00223     void TMult(const Vect<T_> &x, Vect<T_> &y) const;
00224 
00230     void Add(size_t i, size_t j, const T_ &x);
00231 
00234     size_t getColHeight(size_t i) const;
00235 
00241     const T_ operator()(size_t i) const { return _a[i-1]; }
00242 
00247     T_ operator()(size_t i, size_t j) const;
00248 
00253     T_ & operator()(size_t i, size_t j);
00254 
00266     void DiagPrescribe(const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u, int flag=0);
00267 
00270     SkMatrix<T_> & operator=(const SkMatrix<T_> &m);
00271 
00274     SkMatrix<T_> & operator=(const T_ &x);
00275 
00278     SkMatrix<T_> & operator+=(const SkMatrix<T_> &m);
00279 
00282     SkMatrix<T_> & operator+=(const T_ &x);
00283 
00286     SkMatrix<T_> & operator*=(const T_ &x);
00287 
00292     int Factor();
00293 
00300     int Solve(Vect<T_> &b);
00301 
00309     int Solve(const Vect<T_> &b, Vect<T_> &x);
00310 
00313     T_ *getArray() const { return _a; }
00314 
00316     T_ getEntry(size_t i, size_t j) const;
00317 
00318  private:
00319 
00320    mutable ErrorInSkMatrix<T_> _e;
00321    int _dof;
00322 
00323    void set(size_t i, size_t j, const T_ &x);
00324 };
00325 
00326 
00328 //                         I M P L E M E N T A T I O N                       //
00330 
00331 template<class T_>
00332 SkMatrix<T_>::SkMatrix() : _dof(0)
00333 {
00334    _fact = false;
00335 #ifdef _OFELI_DEBUG
00336    std::clog << "An instance of class SkMatrix is constructed.\n";
00337    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00338 #endif
00339 }
00340 
00341 
00342 template<class T_>
00343 SkMatrix<T_>::SkMatrix(const class Mesh &mesh, size_t dof)
00344 {
00345    setMesh(mesh,dof);
00346 }
00347 
00348 
00349 template<class T_>
00350 SkMatrix<T_>::SkMatrix(size_t size) : _dof(0)
00351 {
00352    _zero = 0;
00353    _size = size;
00354    _fact = false;
00355    _ch.resize(_size);
00356    _ch[0] = 0;
00357    for (size_t i=1; i<_size; i++)
00358       _ch[i] = _ch[i-1] + i + 1;
00359    _length = _ch[_size-1]+1;
00360    _a.resize(_length);
00361    _aU.resize(_length);
00362    Clear(_a);
00363    Clear(_aU);
00364    _diag.resize(_size);
00365 #ifdef _OFELI_DEBUG
00366    std::clog << "An instance of class SkMatrix is constructed.\n";
00367    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00368 #endif
00369 }
00370 
00371 
00372 template<class T_>
00373 SkMatrix<T_>::SkMatrix(const Vect<size_t> &ColHt) : _dof(0)
00374 {
00375    _zero = 0;
00376    _ch.resize(_size);
00377    _ch[0] = 0;
00378    for (size_t i=1; i<_size; i++)
00379       _ch[i] = _ch[i-1] + ColHt[i];
00380    _length = _ch[_size-1] + 1;
00381    _a.resize(_length);
00382    _aU.resize(_length);
00383    _diag.resize(_size);
00384    Clear(_a);
00385    Clear(_aU);
00386    _fact = false;
00387 #ifdef _OFELI_DEBUG
00388    std::clog << "An instance of class SkMatrix is constructed.\n";
00389    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00390 #endif
00391 }
00392 
00393 
00394 template<class T_>
00395 SkMatrix<T_>::SkMatrix(const SkMatrix<T_> &m)
00396 {
00397    _size = m._size;
00398    _length = m._length;
00399    _ch.resize(_size);
00400    _ch = m._ch;
00401    _diag.resize(_size);
00402    _diag = m._diag;
00403    _a.resize(_length);
00404    _aU.resize(_length);
00405    _a = m._a;
00406    _aU = m._aU;
00407    _fact = m._fact;
00408    _dof = m._dof;
00409    _zero = T_(0);
00410 #ifdef _OFELI_DEBUG
00411    std::clog << "An instance of class SkMatrix is constructed.\n";
00412    std::clog << "File: " << __FILE__ << ", Line: " << __LINE__ << endl;
00413 #endif
00414 }
00415 
00416 
00417 template<class T_>
00418 SkMatrix<T_>::~SkMatrix()
00419 {
00420 #ifdef _OFELI_DEBUG
00421    std::clog << "An instance of class SkMatrix is destructed.\n";
00422    std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
00423 #endif
00424 }
00425 
00426 
00427 template<class T_>
00428 void SkMatrix<T_>::setMesh(const class Mesh &mesh, size_t dof, size_t type)
00429 {
00430    Matrix<T_>::_init_set_mesh(mesh,dof,type);
00431    if (_dof_type==NODE_DOF)
00432       if (dof)
00433          _length = NodeSkyline(mesh,_ch,dof);
00434       else
00435          _length = NodeSkyline(mesh,_ch);
00436    else if (_dof_type==SIDE_DOF)
00437       if (dof)
00438          _length = SideSkyline(mesh,_ch,dof);
00439       else
00440          _length = SideSkyline(mesh,_ch);
00441    else if (_dof_type==ELEMENT_DOF)
00442       if (dof)
00443          _length = ElementSkyline(mesh,_ch,dof);
00444       else
00445          _length = ElementSkyline(mesh,_ch);
00446    else;
00447       _diag.resize(_size);
00448    _a.resize(_length);
00449    _aU.resize(_length);
00450    _ch[0] = 0;
00451    for (size_t i=1; i<_size; i++)
00452       _ch[i] += _ch[i-1];
00453    Clear(_a);
00454    Clear(_aU);
00455    _zero = T_(0);
00456    _fact = false;
00457 }
00458 
00459 
00460 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00461 template<class T_>
00462 void SkMatrix<T_>::setMesh(size_t dof, const class Mesh &mesh, int code)
00463 {
00464 // This is just to avoid warning on unused variable
00465    dof = 0;
00466    if (mesh.getDim()==0) { }
00467    code = 0;
00468 }
00469 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00470 
00471 
00472 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00473 template<class T_>
00474 void SkMatrix<T_>::setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh)
00475 {
00476 // This is just to avoid warning on unused variable
00477    dof = 0;
00478    nb_eq = 0;
00479    if (mesh.getDim()==0) { }
00480 }
00481 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00482 
00483 
00484 template<class T_>
00485 void SkMatrix<T_>::setSkyline(const class Mesh &mesh)
00486 {
00487    size_t i, k;
00488    _zero = 0;
00489    bool set_nodes = mesh.NodesAreDOF();
00490    bool set_sides = mesh.SidesAreDOF();
00491    bool set_elements = mesh.ElementsAreDOF();
00492    _size = mesh.getNbEq();
00493    if (_dof)
00494       _size = mesh.getNbNodes();
00495    if (set_sides) {
00496       if (_dof) {
00497          _size = mesh.getNbSides();
00498          _length = SideSkyline(mesh,_ch,_dof);
00499       }
00500       _length = SideSkyline(mesh,_ch);
00501    }
00502    else {
00503       if (_dof) {
00504          _size = mesh.getNbNodes();
00505          _length = SideSkyline(mesh,_ch,_dof);
00506       }
00507       _length = SideSkyline(mesh,_ch);
00508    }
00509    _diag.resize(_size);
00510    _ch[0] = 0;
00511    for (i=1; i<_size; i++)
00512       _ch[i] += _ch[i-1];
00513    _a.resize(_length);
00514    _aU.resize(_length);
00515    Clear(_a);
00516    Clear(_aU);
00517    _fact = false;
00518 }
00519 
00520 
00521 template<class T_>
00522 void SkMatrix<T_>::setDiag()
00523 {
00524    for (size_t i=0; i<_size; i++)
00525       _diag[i] = _aU[_ch[i]];
00526 }
00527 
00528 
00529 template<class T_>
00530 void SkMatrix<T_>::Set(size_t i, size_t j, const T_ &x)
00531 {
00532    try {
00533       int k=0, l=0;
00534       if (i>1)
00535          k = int(j-i+_ch[i-1]-_ch[i-2]-1);
00536       if (j>1)
00537          l = int(i-j+_ch[j-1]-_ch[j-2]-1);
00538       if (k>=0 && i>j)
00539          _a[_ch[i-1]+j-i] = x;
00540       else if (l>=0 && i<=j)
00541          _aU[_ch[j-1]+i-j] = x;
00542       else
00543          throw(51);
00544    }
00545    catch(int n) { _e.Message(__FILE__,__LINE__,n,int(i),int(j)); }
00546 }
00547 
00548 
00549 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00550 template<class T_>
00551 void SkMatrix<T_>::SSet(size_t i, size_t j, const T_ &x)
00552 {
00553    int k=0, l=0;
00554    if (i>1)
00555       k = j-i+_ch[i-1]-_ch[i-2]-1;
00556    if (j>1)
00557       l = i-j+_ch[j-1]-_ch[j-2]-1;
00558    if (k>=0 && i>j)
00559       _a[_ch[i-1]+j-i] = x;
00560    else if (l>=0 && i<=j)
00561       _aU[_ch[j-1]+i-j] = x;
00562 }
00563 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00564 
00565 
00566 template<class T_>
00567 void SkMatrix<T_>::MultAdd(const Vect<T_> &x, Vect<T_> &y) const
00568 {
00569    for (size_t i=0; i<_size; i++)
00570       for (size_t j=0; j<_size; j++)
00571          y[i] += operator()(i+1,j+1)*x[j];
00572 }
00573 
00574 
00575 template<class T_>
00576 void SkMatrix<T_>::TMultAdd(const Vect<T_> &x, Vect<T_> &y) const
00577 {
00578    cerr << "TMultAdd is not implemented for class SkMatrix" << endl;
00579 }
00580 
00581 
00582 template<class T_>
00583 void SkMatrix<T_>::MultAdd(T_ a, const Vect<T_> &v, Vect<T_> &w) const
00584 {
00585    for (size_t i=0; i<_size; i++)
00586       for (size_t j=0; j<_size; j++)
00587          w[i] += a * operator()(i+1,j+1)*v[j];
00588 }
00589 
00590 
00591 template<class T_>
00592 void SkMatrix<T_>::Mult(const Vect<T_> &x, Vect<T_> &y) const
00593 {
00594    Clear(y);
00595    MultAdd(x,y);
00596 }
00597 
00598 
00599 template<class T_>
00600 void SkMatrix<T_>::TMult(const Vect<T_> &x, Vect<T_> &y) const
00601 {
00602    Clear(y);
00603    TMultAdd(x,y);
00604 }
00605 
00606 
00607 template<class T_>
00608 void SkMatrix<T_>::Add(size_t i, size_t j, const T_ &x)
00609 {
00610    if (i>j)
00611       _a[_ch[i-1]+j-i] += x;
00612    else if (i<=j)
00613       _aU[_ch[j-1]+i-j] += x;
00614 }
00615 
00616 
00617 template<class T_>
00618 size_t SkMatrix<T_>::getColHeight(size_t i) const
00619 {
00620 #ifdef _BOUNDS_DEBUG
00621    assert(i>0);
00622    assert(i<=_size);
00623 #endif
00624    if (i==1)
00625       return 1;
00626    else
00627       return _ch[i-1]-_ch[i-2];
00628 }
00629 
00630 
00631 template<class T_>
00632 T_ SkMatrix<T_>::operator()(size_t i, size_t j) const
00633 {
00634 #ifdef _BOUNDS_DEBUG
00635    assert(i>0);
00636    assert(i<=_size);
00637    assert(j>0);
00638    assert(j<=_size);
00639 #endif
00640    int k=0, l=0;
00641    if (i>1)
00642       k = int(j-i+_ch[i-1]-_ch[i-2]-1);
00643    if (j>1)
00644       l = int(i-j+_ch[j-1]-_ch[j-2]-1);
00645    if (k>=0 && i>j)
00646       return _a[_ch[i-1]+j-i];
00647    else if (l>=0 && i<=j)
00648       return _aU[_ch[j-1]+i-j];
00649    else
00650       return _zero;
00651 }
00652 
00653 
00654 template<class T_>
00655 T_ & SkMatrix<T_>::operator()(size_t i, size_t j)
00656 {
00657 #ifdef _BOUNDS_DEBUG
00658    assert(i>0);
00659    assert(i<=_size);
00660    assert(j>0);
00661    assert(j<=_size);
00662 #endif
00663    int k=0, l=0;
00664    try {
00665       if (i>1)
00666          k = int(j-i+_ch[i-1]-_ch[i-2]-1);
00667       if (j>1)
00668          l = int(i-j+_ch[j-1]-_ch[j-2]-1);
00669       if (k>=0 && i>j)
00670          return _a[_ch[i-1]+j-i];
00671       else if (l>=0 && i<=j)
00672          return _aU[_ch[j-1]+i-j];
00673       else
00674          throw(81);
00675    }
00676    catch(int n) { _e.Message(__FILE__,__LINE__,n,int(i),int(j)); }
00677    return _temp;
00678 }
00679 
00680 
00681 template<class T_>
00682 void SkMatrix<T_>::DiagPrescribe(const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u, int flag)
00683 {
00684    double p = 0;
00685    for (size_t l=0; l<_size; l++)
00686       p = max(p,_aU[_ch[l]]);
00687    const Node *nd;
00688    for (mesh.topNode(); (nd=mesh.getNode());) {
00689       for (size_t i=1; i<=nd->getNbDOF(); ++i) {
00690          if (nd->getCode(i)>0) {
00691             size_t ii = nd->getDOF(i)-1;
00692             for (size_t j=ii+1+_ch[ii-1]-_ch[ii]; j<=ii; j++) {
00693                b[ii] = p*u[ii];
00694                _a[_ch[ii]+j-ii] = _aU[_ch[ii]+j-ii] = 0;
00695             }
00696             _diag[ii] = _aU[_ch[ii]] = p;
00697          }
00698       }
00699    }
00700 }
00701 
00702 
00703 template<class T_>
00704 SkMatrix<T_> & SkMatrix<T_>::operator=(const SkMatrix<T_> &m)
00705 {
00706    _a = m._a;
00707    _aU = m._aU;
00708    return *this;
00709 }
00710 
00711 
00712 template<class T_>
00713 SkMatrix<T_> & SkMatrix<T_>::operator=(const T_ &x)
00714 {
00715    for (size_t i=0; i<_length; i++)
00716       _a[i] = _aU[i] = x;
00717    return *this;
00718 }
00719 
00720 
00721 template<class T_>
00722 SkMatrix<T_> & SkMatrix<T_>::operator+=(const SkMatrix<T_> &m)
00723 {
00724 #ifdef _BOUNDS_DEBUG
00725    assert(m._length >= _size);
00726 #endif
00727    for (size_t i=0; i<_length; i++) {
00728       _a[i] += m._a[i];
00729       _aU[i] += m._aU[i];
00730    }
00731    return *this;
00732 }
00733 
00734 
00735 template<class T_>
00736 SkMatrix<T_> & SkMatrix<T_>::operator+=(const T_ &x)
00737 {
00738 #ifdef _BOUNDS_DEBUG
00739    assert(m.getSize() >= _size);
00740 #endif
00741    for (size_t i=0; i<_length; i++) {
00742       _a[i] += x;
00743       _aU[i] += x;
00744    }
00745    return *this;
00746 }
00747 
00748 
00749 template<class T_>
00750 SkMatrix<T_> & SkMatrix<T_>::operator*=(const T_ &x)
00751 {
00752    for (size_t i=0; i<_length; i++) {
00753       _a[i] *= x;
00754       _aU[i] *= x;
00755    }
00756    return *this;
00757 }
00758 
00759 
00760 template<class T_>
00761 int SkMatrix<T_>::Factor()
00762 {
00763    size_t k, di, dij, i=0;
00764    if (Abs(_aU[_ch[0]]) == 0)
00765       _e.Message(__FILE__,__LINE__,61,1,0,_aU[_ch[i]]);
00766    for (i=1; i<_size; i++) {
00767       size_t dj = 0;
00768       for (size_t j=di=i+1+_ch[i-1]-_ch[i]; j<i; j++) {
00769          if (j>0)
00770             dj = j+1+_ch[j-1]-_ch[j];
00771          dij = max(di,dj);
00772          for (k=0; k<j-dij; k++)
00773             _a[_ch[i]+j-i] -= _a[_ch[i]+dij+k-i]*_aU[_ch[j]+dij+k-j];
00774          _a[_ch[i]+j-i] /= _aU[_ch[j]];
00775          for (k=0; k<j-dij; k++)
00776             _aU[_ch[i]+j-i] -= _a[_ch[j]+dij+k-j]*_aU[_ch[i]+dij+k-i];
00777       }
00778       for (k=0; k<i-di; k++)
00779          _aU[_ch[i]] -= _a[_ch[i]+k+di-i]*_aU[_ch[i]+k+di-i];
00780       if (Abs(_aU[_ch[i]]) < OFELI_TOLERANCE)
00781          _e.Message(__FILE__,__LINE__,61,int(i+1),0,_aU[_ch[i]]);
00782    }
00783    _fact = true;
00784    return 0;
00785 }
00786 
00787 
00788 template<class T_>
00789 int SkMatrix<T_>::Solve(Vect<T_> &b)
00790 {
00791    size_t di;
00792    try {
00793       if (!_fact)
00794          throw(72);
00795    }
00796    catch(int n) { _e.Message(__FILE__,__LINE__,n); }
00797    size_t i, j;
00798    for (i=1; i<_size; i++) {
00799       di = i+1+_ch[i-1]-_ch[i];
00800       T_ s = 0;
00801       for (j=0; j<i-di; j++)
00802          s += _a[_ch[i]+di+j-i] * b[di+j];
00803       b[i] -= s;
00804    }
00805    for (int k=int(_size-1); k>0; k--) {
00806       if (Abs(_aU[_ch[k]]) < OFELI_TOLERANCE)
00807          _e.Message(__FILE__,__LINE__,71,k+1,0,_aU[_ch[k]]);
00808       b[k] /= _aU[_ch[k]];
00809       di = k+1+_ch[k-1]-_ch[k];
00810       for (j=0; j<k-di; j++)
00811          b[j+di] -= b[k] * _aU[_ch[k]+di+j-k];
00812    }
00813    b[0] /= _aU[_ch[0]];
00814    return 0;
00815 }
00816 
00817 
00818 template<class T_>
00819 int SkMatrix<T_>::Solve(const Vect<T_> &b, Vect<T_> &x)
00820 {
00821    x = b;
00822    return Solve(x);
00823 }
00824 
00825 
00826 template<class T_>
00827 T_ SkMatrix<T_>::getEntry(size_t i, size_t j) const
00828 {
00829    if (i>j)
00830       return _a[_ch[i-1]+j-i];
00831    else
00832       return _aU[_ch[j-1]+i-j];
00833 }
00834 
00835 
00836 template<class T_>
00837 void SkMatrix<T_>::set(size_t i, size_t j, const T_ &x)
00838 {
00839    if (i>j)
00840       return _a[_ch[i-1]+j-i];
00841    else
00842       return _aU[_ch[j-1]+i-j];
00843 }
00844 
00846 //                           ASSOCIATED  FUNCTIONS                           //
00848 
00849 
00853 template<class T_>
00854 ostream& operator<<(ostream& s, const SkMatrix<T_> & a)
00855 {
00856    size_t rmin, rmax, cmin, cmax;
00857    a.getPrintView(rmin,rmax,cmin,cmax);
00858    s.setf(ios::right|ios::scientific);
00859    s << endl;
00860    for (size_t i=rmin; i<=rmax; i++) {
00861       s << "\nRow  " << setw(6) << i << endl;
00862       for (size_t j=cmin; j<=cmax; j++)
00863           s << "  " << setprecision(8) << std::setfill(' ') << setw(18) << a(i,j);
00864       s << endl;
00865    }
00866    return s;
00867 }
00868 
00869 } /* namespace OFELI */
00870 
00871 #endif