NGSolve
4.9
|
00001 #ifndef FILE_PARDISOINVERSE 00002 #define FILE_PARDISOINVERSE 00003 00004 /* *************************************************************************/ 00005 /* File: pardisoinverse.hpp */ 00006 /* Author: Florian Bachinger */ 00007 /* Date: Feb. 04 */ 00008 /* *************************************************************************/ 00009 00010 /* 00011 interface to the PARDISO - package 00012 */ 00013 00014 00015 00017 00018 00019 namespace ngla 00020 { 00021 using ngbla::integer; 00022 00023 00024 /* 00025 template<class TM, 00026 class TV_ROW = typename mat_traits<TM>::TV_ROW, 00027 class TV_COL = typename mat_traits<TM>::TV_COL> 00028 class PardisoInverse : public SparseFactorization 00029 { 00030 integer height; // matrix size in scalars 00031 integer compressed_height; // matrix size after compression in scalars 00032 integer nze, entrysize; 00033 bool print; 00034 00035 integer pt[128]; 00036 integer hparams[64]; 00037 00038 Array<integer,size_t> rowstart, indices; 00039 Array<typename mat_traits<TM>::TSCAL,size_t> matrix; 00040 00041 integer matrixtype; 00042 bool symmetric, spd; 00043 00044 void SetMatrixType(); 00045 00046 bool compressed; 00047 Array<int> compress; 00048 00049 public: 00050 typedef TV_COL TV; 00051 typedef TV_ROW TVX; 00052 typedef typename mat_traits<TM>::TSCAL TSCAL; 00053 00055 PardisoInverse (const SparseMatrix<TM,TV_ROW,TV_COL> & a, 00056 const BitArray * ainner = NULL, 00057 const Array<int> * acluster = NULL, 00058 int symmetric = 0); 00060 00061 template <typename TSUBSET> 00062 void GetPardisoMatrix (const SparseMatrix<TM,TV_ROW,TV_COL> & a, TSUBSET subset); 00063 00064 virtual ~PardisoInverse (); 00066 int VHeight() const { return height/entrysize; } 00068 int VWidth() const { return height/entrysize; } 00070 virtual void Mult (const BaseVector & x, BaseVector & y) const; 00071 00072 00074 virtual ostream & Print (ostream & ost) const; 00075 00076 virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const 00077 { 00078 mu.Append (new MemoryUsageStruct ("SparseChol", nze*sizeof(TM), 1)); 00079 } 00080 virtual BaseVector * CreateVector () const 00081 { 00082 return new VVector<TV> (height/entrysize); 00083 } 00084 }; 00085 */ 00086 00087 00088 00089 00090 00091 00092 00093 00094 00095 template<class TM> 00096 class PardisoInverseTM : public SparseFactorization 00097 { 00098 protected: 00099 integer height; // matrix size in scalars 00100 integer compressed_height; // matrix size after compression in scalars 00101 integer nze, entrysize; 00102 bool print; 00103 00104 integer pt[128]; 00105 integer hparams[64]; 00106 00107 Array<integer,size_t> rowstart, indices; 00108 Array<typename mat_traits<TM>::TSCAL,size_t> matrix; 00109 00110 integer matrixtype; 00111 bool symmetric, spd; 00112 00113 void SetMatrixType(); 00114 00115 bool compressed; 00116 Array<int> compress; 00117 00118 public: 00119 typedef typename mat_traits<TM>::TSCAL TSCAL; 00120 00122 PardisoInverseTM (const SparseMatrixTM<TM> & a, 00123 const BitArray * ainner = NULL, 00124 const Array<int> * acluster = NULL, 00125 int symmetric = 0); 00127 00128 template <typename TSUBSET> 00129 void GetPardisoMatrix (const SparseMatrixTM<TM> & a, TSUBSET subset); 00130 00131 virtual ~PardisoInverseTM (); 00133 int VHeight() const { return height/entrysize; } 00135 int VWidth() const { return height/entrysize; } 00137 virtual ostream & Print (ostream & ost) const; 00138 00139 virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const 00140 { 00141 mu.Append (new MemoryUsageStruct ("SparseChol", nze*sizeof(TM), 1)); 00142 } 00143 }; 00144 00145 00146 template<class TM, 00147 class TV_ROW = typename mat_traits<TM>::TV_ROW, 00148 class TV_COL = typename mat_traits<TM>::TV_COL> 00149 class PardisoInverse : public PardisoInverseTM<TM> 00150 { 00151 using PardisoInverseTM<TM>::height; 00152 using PardisoInverseTM<TM>::compressed_height; 00153 using PardisoInverseTM<TM>::entrysize; 00154 using PardisoInverseTM<TM>::pt; 00155 using PardisoInverseTM<TM>::hparams; 00156 using PardisoInverseTM<TM>::matrixtype; 00157 using PardisoInverseTM<TM>::rowstart; 00158 using PardisoInverseTM<TM>::indices; 00159 using PardisoInverseTM<TM>::matrix; 00160 using PardisoInverseTM<TM>::compressed; 00161 using PardisoInverseTM<TM>::compress; 00162 00163 public: 00164 typedef TV_COL TV; 00165 typedef TV_ROW TVX; 00166 typedef typename mat_traits<TM>::TSCAL TSCAL; 00167 00169 PardisoInverse (const SparseMatrix<TM,TV_ROW,TV_COL> & a, 00170 const BitArray * ainner = NULL, 00171 const Array<int> * acluster = NULL, 00172 int symmetric = 0) 00173 : PardisoInverseTM<TM> (a, ainner, acluster, symmetric) { ; } 00174 00175 virtual ~PardisoInverse () { ; } 00177 virtual void Mult (const BaseVector & x, BaseVector & y) const; 00179 virtual BaseVector * CreateVector () const 00180 { 00181 return new VVector<TV> (height/entrysize); 00182 } 00183 }; 00184 00185 00186 00187 00188 00189 } 00190 00191 #endif