NGSolve  4.9
linalg/sparsematrix.hpp
00001 #ifndef FILE_NGS_SPARSEMATRIX
00002 #define FILE_NGS_SPARSEMATRIX
00003 
00004 /**************************************************************************/
00005 /* File:   sparsematrix.hpp                                               */
00006 /* Author: Joachim Schoeberl                                              */
00007 /* Date:   01. Oct. 94, 15 Jan. 02                                        */
00008 /**************************************************************************/
00009 
00010 namespace ngla
00011 {
00012 
00013   template<class TM>
00014   class SparseMatrixTM ;
00015 
00016   template<class TM, 
00017            class TV_ROW = typename mat_traits<TM>::TV_ROW, 
00018            class TV_COL = typename mat_traits<TM>::TV_COL>
00019   class SparseMatrix;
00020 
00021 
00022   template<class TM, 
00023            class TV = typename mat_traits<TM>::TV_ROW>
00024   class SparseMatrixSymmetric;
00025 
00026 
00027 
00028   class BaseJacobiPrecond;
00029 
00030   template<class TM, 
00031            class TV_ROW = typename mat_traits<TM>::TV_ROW, 
00032            class TV_COL = typename mat_traits<TM>::TV_COL>
00033   class JacobiPrecond;
00034 
00035   template<class TM, 
00036            class TV = typename mat_traits<TM>::TV_ROW>
00037   class JacobiPrecondSymmetric;
00038 
00039 
00040   class BaseBlockJacobiPrecond;
00041 
00042   template<class TM, 
00043            class TV_ROW = typename mat_traits<TM>::TV_ROW, 
00044            class TV_COL = typename mat_traits<TM>::TV_COL>
00045   class BlockJacobiPrecond;
00046 
00047   template<class TM, 
00048            class TV = typename mat_traits<TM>::TV_ROW>
00049   class BlockJacobiPrecondSymmetric;
00050 
00051 
00052 
00053 
00054 
00055   
00056 #ifdef USE_PARDISO
00057   const INVERSETYPE default_inversetype = PARDISO;
00058 #else
00059 #ifdef USE_MUMPS
00060   const INVERSETYPE default_inversetype = MUMPS;
00061 #else
00062   const INVERSETYPE default_inversetype = SPARSECHOLESKY;
00063 #endif
00064 #endif
00065 
00066 
00067 
00068 
00069 
00070 
00074   class NGS_DLL_HEADER MatrixGraph
00075   {
00076   protected:
00078     int size;
00080     int width;
00082     size_t nze; 
00083 
00085     // DynamicMem<int> colnr;
00086     Array<int, size_t> colnr;
00087 
00089     // DynamicMem<size_t> firsti;
00090     Array<size_t, size_t> firsti;
00091   
00093     Array<int> same_nze;
00094 
00096     bool owner;
00097 
00098   public:
00100     MatrixGraph (const Array<int> & elsperrow, int awidth);
00102     MatrixGraph (int as, int max_elsperrow);    
00104     MatrixGraph (const MatrixGraph & graph, bool stealgraph);
00106     MatrixGraph (int size, const Table<int> & rowelements, const Table<int> & colelements, bool symmetric);
00108     // MatrixGraph (const Table<int> & dof2dof, bool symmetric);
00109     virtual ~MatrixGraph ();
00110 
00112     void Compress();
00113   
00115     int GetPosition (int i, int j) const;
00116 
00118     int GetPositionTest (int i, int j) const;
00119 
00121     void GetPositionsSorted (int row, int n, int * pos) const;
00122 
00124     int CreatePosition (int i, int j);
00125 
00126     int Size() const { return size; }
00127 
00128     size_t NZE() const { return nze; }
00129 
00130     FlatArray<int> GetRowIndices(int i) const
00131     // { return FlatArray<int> (int(firsti[i+1]-firsti[i]), colnr+firsti[i]); }
00132     { return FlatArray<int> (int(firsti[i+1]-firsti[i]), &colnr[firsti[i]]); }
00133 
00134     size_t First (int i) const { return firsti[i]; }
00135 
00136     void FindSameNZE();
00137 
00138     ostream & Print (ostream & ost) const;
00139 
00140     virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const;
00141   };
00142 
00143 
00144   
00145 
00146 
00147 
00148 
00149 
00151   class NGS_DLL_HEADER BaseSparseMatrix : virtual public BaseMatrix, 
00152                                           public MatrixGraph
00153   {
00154   protected:
00156     mutable INVERSETYPE inversetype; //  = default_inversetype;    // C++11
00157 
00158   public:
00159     
00160     BaseSparseMatrix (int as, int max_elsperrow)
00161       : MatrixGraph (as, max_elsperrow), inversetype(default_inversetype)
00162     { ; }
00163     
00164     BaseSparseMatrix (const Array<int> & elsperrow, int awidth)
00165       : MatrixGraph (elsperrow, awidth), inversetype(default_inversetype)
00166     { ; }
00167 
00168     BaseSparseMatrix (int size, const Table<int> & rowelements, 
00169                       const Table<int> & colelements, bool symmetric)
00170       : MatrixGraph (size, rowelements, colelements, symmetric),
00171         inversetype(default_inversetype)
00172     { ; }
00173 
00174     BaseSparseMatrix (const MatrixGraph & agraph, bool stealgraph)
00175       : MatrixGraph (agraph, stealgraph), inversetype(default_inversetype)
00176     { ; }   
00177 
00178     BaseSparseMatrix (const BaseSparseMatrix & amat)
00179       : BaseMatrix(amat), MatrixGraph (amat, 0), inversetype(default_inversetype)
00180     { ; }   
00181 
00182 
00183     virtual ~BaseSparseMatrix ();
00184 
00185     BaseSparseMatrix & operator= (double s)
00186     {
00187       AsVector() = s;
00188       return *this;
00189     }
00190 
00191     BaseSparseMatrix & Add (double s, const BaseSparseMatrix & m2)
00192     {
00193       AsVector() += s * m2.AsVector();
00194       return *this;
00195     }
00196 
00197     virtual BaseJacobiPrecond * CreateJacobiPrecond (const BitArray * inner = 0) const 
00198     {
00199       throw Exception ("BaseSparseMatrix::CreateJacobiPrecond");
00200     }
00201 
00202     virtual BaseBlockJacobiPrecond * 
00203     CreateBlockJacobiPrecond (Table<int> & blocks,
00204                               const BaseVector * constraint = 0,
00205                               const ngcomp::Preconditioner * acoarsegridprecond = 0, 
00206                               bool parallel  = 1,
00207                               const BitArray * freedofs = NULL) const
00208     { 
00209       throw Exception ("BaseSparseMatrix::CreateBlockJacobiPrecond");
00210     }
00211 
00212 
00213     virtual BaseMatrix * 
00214     InverseMatrix (const BitArray * subset = 0) const
00215     { 
00216       throw Exception ("BaseSparseMatrix::CreateInverse called");
00217     }
00218 
00219     virtual BaseMatrix * 
00220     InverseMatrix (const Array<int> * clusters) const
00221     { 
00222       throw Exception ("BaseSparseMatrix::CreateInverse called");
00223     }
00224 
00225     virtual BaseSparseMatrix * Restrict (const SparseMatrixTM<double> & prol,
00226                                          BaseSparseMatrix* cmat = NULL ) const
00227     {
00228       throw Exception ("BaseSparseMatrix::Restrict");
00229     }
00230 
00231     //virtual void AllReduceHO () const = 0;
00232 
00233 
00234     virtual INVERSETYPE SetInverseType ( INVERSETYPE ainversetype ) const
00235     {
00236 
00237       INVERSETYPE old_invtype = inversetype;
00238       inversetype = ainversetype; 
00239       return old_invtype;
00240     }
00241 
00242     virtual INVERSETYPE SetInverseType ( string ainversetype ) const;
00243 
00244     virtual INVERSETYPE  GetInverseType () const
00245     { return inversetype; }
00246 
00247 
00248 
00249   };
00250 
00252   template<class TM>
00253   class  NGS_DLL_HEADER SparseMatrixTM : public BaseSparseMatrix, 
00254                                          public S_BaseMatrix<typename mat_traits<TM>::TSCAL>
00255   {
00256   protected:
00257     // DynamicMem<TM> data;
00258     Array<TM, size_t> data;
00259     VFlatVector<typename mat_traits<TM>::TSCAL> asvec;
00260     TM nul;
00261 
00262   public:
00263     typedef typename mat_traits<TM>::TSCAL TSCAL;
00264 
00265     SparseMatrixTM (int as, int max_elsperrow);
00266     SparseMatrixTM (const Array<int> & elsperrow, int awidth);
00267     SparseMatrixTM (int size, const Table<int> & rowelements, 
00268                     const Table<int> & colelements, bool symmetric);
00269     SparseMatrixTM (const MatrixGraph & agraph, bool stealgraph);
00270     SparseMatrixTM (const SparseMatrixTM & amat);
00271     virtual ~SparseMatrixTM ();
00272 
00273     int Height() const { return size; }
00274     int Width() const { return width; }
00275     virtual int VHeight() const { return size; }
00276     virtual int VWidth() const { return width; }
00277 
00278     TM & operator[] (int i)  { return data[i]; }
00279     const TM & operator[] (int i) const { return data[i]; }
00280 
00281     TM & operator() (int row, int col)
00282     {
00283       return data[CreatePosition(row, col)];
00284     }
00285 
00286     const TM & operator() (int row, int col) const
00287     {
00288       int pos = GetPositionTest (row,col);
00289       if (pos != -1)
00290         return data[pos];
00291       else
00292         return nul;
00293     }
00294 
00295     FlatVector<TM> GetRowValues(int i) const
00296     { return FlatVector<TM> (firsti[i+1]-firsti[i], &data[firsti[i]]); }
00297 
00298 
00299     virtual void AddElementMatrix(const FlatArray<int> & dnums1, 
00300                                   const FlatArray<int> & dnums2, 
00301                                   const FlatMatrix<TSCAL> & elmat);
00302 
00303     virtual BaseVector & AsVector() 
00304     {
00305       asvec.AssignMemory (nze*sizeof(TM)/sizeof(TSCAL), (void*)&data[0]);
00306       return asvec; 
00307     }
00308 
00309     virtual const BaseVector & AsVector() const
00310     { 
00311       const_cast<VFlatVector<TSCAL>&> (asvec).
00312         AssignMemory (nze*sizeof(TM)/sizeof(TSCAL), (void*)&data[0]);
00313       return asvec; 
00314     }
00315 
00317     virtual ostream & Print (ostream & ost) const;
00318 
00320     virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const;
00321 
00322 
00323   };
00324   
00325 
00326 
00327 
00328   template<class TM, class TV_ROW, class TV_COL>
00329   class NGS_DLL_HEADER SparseMatrix : virtual public SparseMatrixTM<TM>
00330   {
00331   public:
00332     using SparseMatrixTM<TM>::firsti;
00333     using SparseMatrixTM<TM>::colnr;
00334     using SparseMatrixTM<TM>::data;
00335 
00336     // using MatrixGraph::firsti;
00337 
00338 
00339     typedef typename mat_traits<TM>::TSCAL TSCAL;
00340     typedef TV_ROW TVX;
00341     typedef TV_COL TVY;
00342     typedef typename mat_scale_type<TVX,Complex>::TMAT TVXC;
00343     typedef typename mat_scale_type<TVY,Complex>::TMAT TVYC;
00344 
00346     SparseMatrix (int as, int max_elsperrow)
00347       : SparseMatrixTM<TM> (as, max_elsperrow) { ; }
00348 
00349     SparseMatrix (const Array<int> & aelsperrow)
00350       : SparseMatrixTM<TM> (aelsperrow, aelsperrow.Size()) { ; }
00351 
00352     SparseMatrix (const Array<int> & aelsperrow, int awidth)
00353       : SparseMatrixTM<TM> (aelsperrow, awidth) { ; }
00354 
00355     SparseMatrix (int size, const Table<int> & rowelements, 
00356                   const Table<int> & colelements, bool symmetric)
00357       : SparseMatrixTM<TM> (size, rowelements, colelements, symmetric) { ; }
00358 
00359     SparseMatrix (const MatrixGraph & agraph, bool stealgraph)
00360       : SparseMatrixTM<TM> (agraph, stealgraph) { ; }
00361 
00362     SparseMatrix (const SparseMatrix & amat)
00363       : SparseMatrixTM<TM> (amat) { ; }
00364 
00365     SparseMatrix (const SparseMatrixTM<TM> & amat)
00366       : SparseMatrixTM<TM> (amat) { ; }
00367 
00368     virtual BaseMatrix * CreateMatrix () const;
00369     // virtual BaseMatrix * CreateMatrix (const Array<int> & elsperrow) const;
00371     virtual BaseVector * CreateVector () const;
00372 
00373     virtual BaseJacobiPrecond * 
00374     CreateJacobiPrecond (const BitArray * inner) const
00375     { 
00376       return new JacobiPrecond<TM,TV_ROW,TV_COL> (*this, inner);
00377     }
00378 
00379     virtual BaseBlockJacobiPrecond * 
00380     CreateBlockJacobiPrecond (Table<int> & blocks,
00381                               const BaseVector * constraint = 0,
00382                               const ngcomp::Preconditioner * acoarsegridprecond = 0, 
00383                               bool parallel  = 1,
00384                               const BitArray * freedofs = NULL) const
00385     { 
00386       return new BlockJacobiPrecond<TM,TV_ROW,TV_COL> (*this, blocks );
00387     }
00388 
00389     virtual BaseMatrix * InverseMatrix (const BitArray * subset = 0) const;
00390     virtual BaseMatrix * InverseMatrix (const Array<int> * clusters) const;
00391 
00392     virtual BaseSparseMatrix * Restrict (const SparseMatrixTM<double> & prol,
00393                                          BaseSparseMatrix* cmat = NULL ) const
00394     { return 0; }
00395 
00396 
00397   
00399     // template<class TEL>
00400     inline TVY RowTimesVector (int row, const FlatVector<TVX> & vec) const
00401     {
00402       typedef typename mat_traits<TVY>::TSCAL TTSCAL;
00403       TVY sum = TTSCAL(0);
00404       for (size_t j = firsti[row]; j < firsti[row+1]; j++)
00405         sum += data[j] * vec(colnr[j]);
00406       return sum;
00407 
00408       /*
00409       int nj = firsti[row+1] - firsti[row];
00410       const int * colpi = colnr+firsti[row];
00411       const TM * datap = data+firsti[row];
00412       const TVX * vecp = vec.Addr(0);
00413     
00414       typedef typename mat_traits<TVY>::TSCAL TTSCAL;
00415       TVY sum = TTSCAL(0);
00416       for (int j = 0; j < nj; j++, colpi++, datap++)
00417         sum += *datap * vecp[*colpi];
00418       return sum;
00419       */
00420     }
00421 
00423     void AddRowTransToVector (int row, TVY el, FlatVector<TVX> & vec) const
00424     {
00425       // cannot be optimized by the compiler, since hvec could alias
00426       // the matrix --> keyword 'restrict' will hopefully solve this problem
00427       size_t first = firsti[row];
00428       size_t last = firsti[row+1];
00429       for (size_t j = first; j < last; j++)
00430         vec(colnr[j]) += Trans(data[j]) * el;
00431       /*
00432       int nj = firsti[row+1] - firsti[row];
00433       const int * colpi = colnr+firsti[row];
00434       const TM * datap = data+firsti[row];
00435       TVX * vecp = vec.Addr(0);
00436       for (int j = 0; j < nj; j++, colpi++, datap++)
00437         vecp[*colpi] += Trans(*datap) * el; 
00438       */
00439     }
00440 
00441 
00442     virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const;
00443     virtual void MultTransAdd (double s, const BaseVector & x, BaseVector & y) const;
00444     virtual void MultAdd (Complex s, const BaseVector & x, BaseVector & y) const;
00445     virtual void MultTransAdd (Complex s, const BaseVector & x, BaseVector & y) const;
00446   };
00447 
00448 
00450   template<class TM>
00451   class NGS_DLL_HEADER SparseMatrixSymmetricTM : virtual public SparseMatrixTM<TM>
00452   {
00453   protected:
00454     SparseMatrixSymmetricTM (int as, int max_elsperrow)
00455       : SparseMatrixTM<TM> (as, max_elsperrow) { ; }
00456 
00457     SparseMatrixSymmetricTM (const Array<int> & elsperrow)
00458       : SparseMatrixTM<TM> (elsperrow, elsperrow.Size()) { ; }
00459 
00460     SparseMatrixSymmetricTM (int size, const Table<int> & rowelements)
00461       : SparseMatrixTM<TM> (size, rowelements, rowelements, true) { ; }
00462 
00463     SparseMatrixSymmetricTM (const MatrixGraph & agraph, bool stealgraph)
00464       : SparseMatrixTM<TM> (agraph, stealgraph) { ; }
00465 
00466     SparseMatrixSymmetricTM (const SparseMatrixSymmetricTM & amat)
00467       : SparseMatrixTM<TM> (amat) { ; }
00468 
00469   public:
00470     typedef typename mat_traits<TM>::TSCAL TSCAL;
00471     virtual void AddElementMatrix(const FlatArray<int> & dnums, const FlatMatrix<TSCAL> & elmat);
00472 
00473     virtual void AddElementMatrix(const FlatArray<int> & dnums1, 
00474                                   const FlatArray<int> & dnums2, 
00475                                   const FlatMatrix<TSCAL> & elmat)
00476     {
00477       AddElementMatrix (dnums1, elmat);
00478     }
00479   };
00480 
00481 
00483   template<class TM, class TV>
00484   class NGS_DLL_HEADER SparseMatrixSymmetric : virtual public SparseMatrixSymmetricTM<TM>, 
00485                                                virtual public SparseMatrix<TM,TV,TV>
00486   {
00487 
00488   public:
00489     using SparseMatrixTM<TM>::firsti;
00490     using SparseMatrixTM<TM>::colnr;
00491     using SparseMatrixTM<TM>::data;
00492 
00493     typedef typename mat_traits<TM>::TSCAL TSCAL;
00494     typedef TV TV_COL;
00495     typedef TV TV_ROW;
00496     typedef TV TVY;
00497     typedef TV TVX;
00498 
00499     SparseMatrixSymmetric (int as, int max_elsperrow)
00500       : SparseMatrixTM<TM> (as, max_elsperrow) , 
00501         SparseMatrixSymmetricTM<TM> (as, max_elsperrow),
00502         SparseMatrix<TM,TV,TV> (as, max_elsperrow)
00503     { ; }
00504   
00505     SparseMatrixSymmetric (const Array<int> & elsperrow)
00506       : SparseMatrixTM<TM> (elsperrow, elsperrow.Size()), 
00507         SparseMatrixSymmetricTM<TM> (elsperrow),
00508         SparseMatrix<TM,TV,TV> (elsperrow, elsperrow.Size())
00509     { ; }
00510 
00511     SparseMatrixSymmetric (int size, const Table<int> & rowelements)
00512       : SparseMatrixTM<TM> (size, rowelements, rowelements, true),
00513         SparseMatrixSymmetricTM<TM> (size, rowelements),
00514         SparseMatrix<TM,TV,TV> (size, rowelements, rowelements, true)
00515     { ; }
00516 
00517     SparseMatrixSymmetric (const MatrixGraph & agraph, bool stealgraph)
00518       : SparseMatrixTM<TM> (agraph, stealgraph), 
00519         SparseMatrixSymmetricTM<TM> (agraph, stealgraph),
00520         SparseMatrix<TM,TV,TV> (agraph, stealgraph)
00521     { ; }
00522 
00523 
00524     SparseMatrixSymmetric (const SparseMatrixSymmetric & amat)
00525       : SparseMatrixTM<TM> (amat), 
00526         SparseMatrixSymmetricTM<TM> (amat),
00527         SparseMatrix<TM,TV,TV> (amat)
00528     { 
00529       this->AsVector() = amat.AsVector(); 
00530     }
00531 
00532     SparseMatrixSymmetric (const SparseMatrixSymmetricTM<TM> & amat)
00533       : SparseMatrixTM<TM> (amat), 
00534         SparseMatrixSymmetricTM<TM> (amat),
00535         SparseMatrix<TM,TV,TV> (amat)
00536     { 
00537       this->AsVector() = amat.AsVector(); 
00538     }
00539 
00540   
00541 
00542 
00543 
00544 
00545 
00547     virtual ~SparseMatrixSymmetric ();
00548 
00549     SparseMatrixSymmetric & operator= (double s)
00550     {
00551       this->AsVector() = s;
00552       return *this;
00553     }
00554 
00555     virtual BaseMatrix * CreateMatrix () const
00556     {
00557       return new SparseMatrixSymmetric (*this);
00558     }
00559 
00560     /*
00561     virtual BaseMatrix * CreateMatrix (const Array<int> & elsperrow) const
00562     {
00563       return new SparseMatrix<TM,TV,TV>(elsperrow);
00564     }
00565     */
00566 
00567     virtual BaseJacobiPrecond * CreateJacobiPrecond (const BitArray * inner) const
00568     { 
00569       return new JacobiPrecondSymmetric<TM,TV> (*this, inner);
00570     }
00571 
00572     virtual BaseBlockJacobiPrecond * 
00573     CreateBlockJacobiPrecond (Table<int> & blocks,
00574                               const BaseVector * constraint = 0,
00575                               const ngcomp::Preconditioner * acoarsegridprecond = 0, 
00576                               bool parallel  = 1,
00577                               const BitArray * freedofs = NULL) const
00578     { 
00579       return new BlockJacobiPrecondSymmetric<TM,TV> (*this, blocks);
00580     }
00581 
00582 
00583     virtual BaseSparseMatrix * Restrict (const SparseMatrixTM<double> & prol,
00584                                          BaseSparseMatrix* cmat = NULL ) const;
00585 
00587     virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const;
00588 
00589     virtual void MultTransAdd (double s, const BaseVector & x, BaseVector & y) const
00590     {
00591       MultAdd (s, x, y);
00592     }
00593 
00594 
00595     /*
00596       y += s L * x
00597     */
00598     virtual void MultAdd1 (double s, const BaseVector & x, BaseVector & y,
00599                            const BitArray * ainner = NULL,
00600                            const Array<int> * acluster = NULL) const;
00601 
00602 
00603     /*
00604       y += s (D + L^T) * x
00605     */
00606     virtual void MultAdd2 (double s, const BaseVector & x, BaseVector & y,
00607                            const BitArray * ainner = NULL,
00608                            const Array<int> * acluster = NULL) const;
00609     
00610 
00611 
00612 
00613 
00614     using SparseMatrix<TM,TV,TV>::RowTimesVector;
00615     using SparseMatrix<TM,TV,TV>::AddRowTransToVector;
00616 
00617 
00618     TV_COL RowTimesVectorNoDiag (int row, const FlatVector<TVX> & vec) const
00619     {
00620       size_t last = firsti[row+1];
00621       size_t first = firsti[row];
00622       if (last == first) return TVY(0);
00623       if (colnr[last-1] == row) last--;
00624 
00625       typedef typename mat_traits<TVY>::TSCAL TTSCAL;
00626       TVY sum = TTSCAL(0);
00627 
00628       for (size_t j = first; j < last; j++)
00629         sum += data[j] * vec(colnr[j]);
00630       return sum;
00631 
00632       /*
00633       const int * colpi = colnr+firsti[row];
00634       const TM * datap = data+firsti[row];
00635       const TVX * vecp = vec.Addr(0);
00636 
00637       for (int j = first; j < last; j++, colpi++, datap++)
00638         sum += *datap * vecp[*colpi];
00639       return sum;
00640       */
00641     }
00642 
00643     void AddRowTransToVectorNoDiag (int row, TVY el, FlatVector<TVX> & vec) const
00644     {
00645       size_t first = firsti[row];
00646       size_t last = firsti[row+1];
00647       if (first == last) return;
00648       if (this->colnr[last-1] == row) last--;
00649 
00650       for (size_t j = first; j < last; j++)
00651         vec(colnr[j]) += Trans(data[j]) * el;
00652 
00653       /*
00654         // hand tuning, no difference anymore
00655       size_t last = firsti[row+1];
00656       size_t first = firsti[row];
00657       if (first == last) return;
00658 
00659       if (colnr[last-1] == row) last--;
00660 
00661       int nj = last - firsti[row];
00662       const int * colpi = colnr+firsti[row];
00663       const TM * datap = data+firsti[row];
00664       TVX * vecp = vec.Addr(0);
00665       for (int j = 0; j < nj; j++, colpi++, datap++)
00666         vecp[*colpi] += Trans(*datap) * el;
00667       */
00668     }
00669   
00670     BaseSparseMatrix & AddMerge (double s, const SparseMatrixSymmetric  & m2);
00671 
00672     virtual BaseMatrix * InverseMatrix (const BitArray * subset = 0) const;
00673     virtual BaseMatrix * InverseMatrix (const Array<int> * clusters) const;
00674   };
00675 
00676 
00677 
00678 }
00679 
00680 #endif