NGSolve
4.9
|
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