NGSolve
4.9
|
00001 #ifndef FILE_SPARSECHOLESKY 00002 #define FILE_SPARSECHOLESKY 00003 00004 /* *************************************************************************/ 00005 /* File: sparsecholesky.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 18. Jun. 97 */ 00008 /* *************************************************************************/ 00009 00010 /* 00011 sparse cholesky factorization 00012 */ 00013 00014 namespace ngla 00015 { 00016 00017 class SparseFactorization : public BaseMatrix 00018 { 00019 protected: 00020 const BaseSparseMatrix & matrix; 00021 const BitArray * inner; 00022 const Array<int> * cluster; 00023 bool smooth_is_projection; 00024 00025 public: 00026 SparseFactorization (const BaseSparseMatrix & amatrix, 00027 const BitArray * ainner, 00028 const Array<int> * acluster); 00029 00030 virtual void Smooth (BaseVector & u, const BaseVector & f, BaseVector & y) const; 00031 00032 bool SmoothIsProjection () const { return smooth_is_projection; } 00033 }; 00034 00035 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 SparseCholesky : public SparseFactorization 00046 { 00047 int height, nze; 00048 00049 /* 00050 DynamicMem<int> order, firstinrow, firstinrow_ri, rowindex2, blocknrs; 00051 DynamicMem<TM> lfact; 00052 DynamicMem<TM> diag; 00053 */ 00054 Array<int, size_t> order, firstinrow, firstinrow_ri, rowindex2, blocknrs; 00055 Array<TM, size_t> lfact; 00056 Array<TM, size_t> diag; 00057 00059 MinimumDegreeOrdering * mdo; 00060 int maxrow; 00061 const SparseMatrix<TM,TV_ROW,TV_COL> & mat; 00062 00063 public: 00064 typedef TV_COL TV; 00065 typedef TV_ROW TVX; 00066 typedef typename mat_traits<TV_ROW>::TSCAL TSCAL_VEC; 00067 typedef typename mat_traits<TM>::TSCAL TSCAL_MAT; 00068 00070 SparseCholesky (const SparseMatrix<TM,TV_ROW,TV_COL> & a, 00071 const BitArray * ainner = NULL, 00072 const Array<int> * acluster = NULL, 00073 bool allow_refactor = 0); 00075 ~SparseCholesky (); 00077 int VHeight() const { return height; } 00079 int VWidth() const { return height; } 00081 void Allocate (const Array<int> & aorder, 00082 const Array<MDOVertex> & vertices, 00083 const int * blocknr); 00085 void Factor (); // const int * blocknr); 00087 void FactorNew (const SparseMatrix<TM,TV_ROW,TV_COL> & a); 00089 virtual void Mult (const BaseVector & x, BaseVector & y) const; 00091 virtual void MultAdd (TSCAL_VEC s, const BaseVector & x, BaseVector & y) const; 00099 // virtual void Smooth (BaseVector & u, const BaseVector & f, BaseVector & y) const; 00101 virtual ostream & Print (ostream & ost) const; 00102 00103 virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const 00104 { 00105 mu.Append (new MemoryUsageStruct ("SparseChol", nze*sizeof(TM), 1)); 00106 } 00107 00108 00110 void Set (int i, int j, const TM & val); 00112 const TM & Get (int i, int j) const; 00114 void SetOrig (int i, int j, const TM & val) 00115 { Set (order[i], order[j], val); } 00116 00117 virtual BaseVector * CreateVector () const 00118 { 00119 return new VVector<TV> (height); 00120 } 00121 }; 00122 00123 } 00124 00125 #endif