NGSolve
4.9
|
00001 #ifndef FILE_BLOCK_JACOBI 00002 #define FILE_BLOCK_JACOBI 00003 00004 /* *************************************************************************/ 00005 /* File: blockjacobi.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 06. Oct. 96 */ 00008 /* *************************************************************************/ 00009 00010 namespace ngla 00011 { 00012 00016 class BaseBlockJacobiPrecond : virtual public BaseMatrix 00017 { 00018 public: 00019 // enum COARSE_TYPE { NO_COARSE = 0, USER_COARSE = 1, DIRECT_COARSE = 2, SMOOTHING_COARSE = 3 }; 00020 protected: 00022 Table<int> & blocktable; 00024 int maxbs; 00025 00026 // COARSE_TYPE ct; 00027 public: 00029 BaseBlockJacobiPrecond (Table<int> & ablocktable); 00030 00032 virtual ~BaseBlockJacobiPrecond (); 00033 00035 virtual void GSSmooth (BaseVector & x, const BaseVector & b, 00036 int steps = 1) const = 0; 00037 00039 virtual void GSSmoothPartial (BaseVector & x, const BaseVector & b, BaseVector & y) const 00040 { 00041 GSSmooth (x, b, 1); 00042 } 00043 00044 00046 virtual void GSSmoothResiduum (BaseVector & x, const BaseVector & b, 00047 BaseVector & res, int steps = 1) const = 0; 00048 00050 virtual void GSSmoothBack (BaseVector & x, const BaseVector & b, 00051 int steps = 1) const = 0; 00052 00053 virtual void GSSmoothBackPartial (BaseVector & x, const BaseVector & b, BaseVector & y) const 00054 { 00055 GSSmoothBack (x, b, 1); 00056 } 00057 00058 00060 int Reorder (FlatArray<int> block, const MatrixGraph & graph, 00061 FlatArray<int> usedflags, // in and out: array of -1, size = graph.size 00062 LocalHeap & lh); 00063 00064 /* 00065 virtual void SetCoarseType ( string act) 00066 { 00067 if ( strcmp ( act.c_str(), "DIRECT_COARSE") == 0 ) 00068 { 00069 ct = DIRECT_COARSE; 00070 } 00071 else if ( strcmp ( act.c_str(), "NO_COARSE") == 0 ) 00072 { 00073 ct = NO_COARSE; 00074 } 00075 else if ( strcmp ( act.c_str(), "USER_COARSE") == 0 ) 00076 { 00077 ct = USER_COARSE; 00078 } 00079 else if ( strcmp ( act.c_str(), "SMOOTHING_COARSE") == 0 ) 00080 { 00081 ct = SMOOTHING_COARSE; 00082 } 00083 else 00084 { 00085 ct = NO_COARSE; 00086 } 00087 } 00088 */ 00089 00090 /* 00091 virtual void InitCoarseType (string act, const BitArray * freedofs) 00092 { 00093 cerr << "BaseBlockJacobiPrecond :: InitCoarseType not implemented!" << endl; 00094 } 00095 */ 00096 }; 00097 00098 00099 00100 00105 template <class TM, class TV_ROW, class TV_COL> 00106 class BlockJacobiPrecond : virtual public BaseBlockJacobiPrecond, 00107 virtual public S_BaseMatrix<typename mat_traits<TM>::TSCAL> 00108 { 00109 protected: 00111 const SparseMatrix<TM,TV_ROW,TV_COL> & mat; 00113 Array<Matrix<TM>*> invdiag; 00114 00115 00116 public: 00117 // typedef typename mat_traits<TM>::TV_ROW TVX; 00118 typedef TV_ROW TVX; 00119 typedef typename mat_traits<TM>::TSCAL TSCAL; 00120 00122 BlockJacobiPrecond (const SparseMatrix<TM,TV_ROW,TV_COL> & amat, 00123 Table<int> & ablocktable); 00125 virtual ~BlockJacobiPrecond (); 00126 00127 int Height() const { return mat.Height(); } 00128 virtual int VHeight() const { return mat.Height(); } 00129 int Width() const { return mat.Width(); } 00130 virtual int VWidth() const { return mat.Width(); } 00131 00133 virtual BaseVector * CreateVector () const 00134 { 00135 return mat.CreateVector(); 00136 } 00137 00138 00140 virtual void MultAdd (TSCAL s, const BaseVector & x, BaseVector & y) const 00141 { 00142 static Timer timer("BlockJacobi::MultAdd"); 00143 RegionTimer reg (timer); 00144 00145 FlatVector<TVX> fx = x.FV<TVX> (); 00146 FlatVector<TVX> fy = y.FV<TVX> (); 00147 00148 Vector<TVX> hxmax(maxbs); 00149 00150 for (int i = 0; i < blocktable.Size(); i++) 00151 { 00152 FlatArray<int> ind = blocktable[i]; 00153 if (!ind.Size()) continue; 00154 00155 FlatVector<TVX> hx(ind.Size(), hxmax.Addr(0)); 00156 00157 hx = s * fx(ind); 00158 fy(ind) += *invdiag[i] * hx; 00159 } 00160 } 00161 00162 00164 virtual void MultTransAdd (TSCAL s, const BaseVector & x, BaseVector & y) const 00165 { 00166 static int timer = NgProfiler::CreateTimer ("BlockJacobi::MultTransAdd"); 00167 NgProfiler::RegionTimer reg (timer); 00168 00169 int i, j; 00170 FlatVector<TVX> fx = x.FV<TVX> (); 00171 FlatVector<TVX> fy = y.FV<TVX> (); 00172 00173 Vector<TVX> hxmax(maxbs); 00174 Vector<TVX> hymax(maxbs); 00175 00176 for (i = 0; i < blocktable.Size(); i++) 00177 { 00178 int bs = blocktable[i].Size(); 00179 if (!bs) continue; 00180 00181 FlatVector<TVX> hx(bs, hxmax.Addr(0)); 00182 FlatVector<TVX> hy(bs, hymax.Addr(0)); 00183 00184 for (j = 0; j < bs; j++) 00185 hx(j) = fx(blocktable[i][j]); 00186 00187 hy = Trans(*invdiag[i]) * hx; 00188 00189 for (j = 0; j < bs; j++) 00190 fy(blocktable[i][j]) += s * hy(j); 00191 } 00192 } 00193 00194 00196 virtual void GSSmooth (BaseVector & x, const BaseVector & b, 00197 int steps = 1) const 00198 { 00199 FlatVector<TVX> fb = b.FV<TVX> (); 00200 FlatVector<TVX> fx = x.FV<TVX> (); 00201 00202 Vector<TVX> hxmax(maxbs); 00203 Vector<TVX> hymax(maxbs); 00204 for (int k = 0; k < steps; k++) 00205 for (int i = 0; i < blocktable.Size(); i++) 00206 { 00207 int bs = blocktable[i].Size(); 00208 if (!bs) continue; 00209 00210 FlatVector<TVX> hx(bs, hxmax.Addr(0)); 00211 FlatVector<TVX> hy(bs, hymax.Addr(0)); 00212 00213 for (int j = 0; j < bs; j++) 00214 { 00215 int jj = blocktable[i][j]; 00216 hx(j) = fb(jj) - mat.RowTimesVector (jj, fx); 00217 } 00218 00219 hy = (*invdiag[i]) * hx; 00220 00221 for (int j = 0; j < bs; j++) 00222 fx(blocktable[i][j]) += hy(j); 00223 } 00224 } 00225 00226 00227 virtual void GSSmoothResiduum (BaseVector & x, const BaseVector & b, 00228 BaseVector & res, int steps = 1) const 00229 { 00230 GSSmooth (x, b, 1); 00231 res = b - mat * x; 00232 } 00233 00234 00235 00236 00238 virtual void GSSmoothBack (BaseVector & x, const BaseVector & b, 00239 int steps = 1) const 00240 { 00241 const FlatVector<TVX> fb = b.FV<TVX> (); 00242 FlatVector<TVX> fx = x.FV<TVX> (); 00243 00244 Vector<TVX> hxmax(maxbs); 00245 Vector<TVX> hymax(maxbs); 00246 00247 for (int k = 0; k < steps; k++) 00248 for (int i = blocktable.Size()-1; i >= 0; i--) 00249 { 00250 int bs = blocktable[i].Size(); 00251 if (!bs) continue; 00252 00253 FlatVector<TVX> hx(bs, hxmax.Addr(0)); 00254 FlatVector<TVX> hy(bs, hymax.Addr(0)); 00255 00256 for (int j = 0; j < bs; j++) 00257 { 00258 int jj = blocktable[i][j]; 00259 hx(j) = fb(jj) - mat.RowTimesVector (jj, fx); 00260 } 00261 00262 hy = (*invdiag[i]) * hx; 00263 00264 for (int j = 0; j < bs; j++) 00265 fx(blocktable[i][j]) += hy(j); 00266 } 00267 } 00268 00270 virtual void GSSmoothNumbering (BaseVector & x, const BaseVector & b, 00271 const Array<int> & numbering, 00272 int forward = 1) const 00273 { 00274 ; 00275 } 00276 00277 virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const 00278 { 00279 int nels = 0; 00280 for (int i = 0; i < blocktable.Size(); i++) 00281 { 00282 int bs = blocktable[i].Size(); 00283 nels += bs*bs; 00284 } 00285 mu.Append (new MemoryUsageStruct ("BlockJac", nels*sizeof(TM), blocktable.Size())); 00286 } 00287 00288 00289 }; 00290 00291 00292 00293 00294 00295 /* **************** SYMMETRIC ****************** */ 00296 00297 00299 template <class TM, class TV> 00300 class BlockJacobiPrecondSymmetric : 00301 virtual public BaseBlockJacobiPrecond, 00302 virtual public S_BaseMatrix<typename mat_traits<TM>::TSCAL> 00303 { 00304 protected: 00305 const SparseMatrixSymmetric<TM,TV> & mat; 00306 00307 // 00308 enum { NBLOCKS = 20 }; 00309 // DynamicMem<int> blockstart, blocksize, blockbw; 00310 // DynamicMem<TM> data[NBLOCKS]; 00311 00312 Array<int, size_t> blockstart, blocksize, blockbw; 00313 Array<TM, size_t> data[NBLOCKS]; 00314 00315 00316 bool lowmem; 00317 public: 00318 00319 typedef TV TVX; 00320 typedef typename mat_traits<TM>::TSCAL TSCAL; 00321 00323 BlockJacobiPrecondSymmetric (const SparseMatrixSymmetric<TM,TV> & amat, 00324 Table<int> & ablocktable); 00325 00326 /* 00327 BlockJacobiPrecondSymmetric (const SparseMatrixSymmetric<TM,TV> & amat, 00328 const FlatVector<TVX> & constraint, 00329 Table<int> & ablocktable); 00330 */ 00332 virtual ~BlockJacobiPrecondSymmetric (); 00333 00334 FlatBandCholeskyFactors<TM> InvDiag (int i) const 00335 { 00336 return FlatBandCholeskyFactors<TM> (blocksize[i], 00337 blockbw[i], 00338 const_cast<TM*>(&data[i%NBLOCKS][blockstart[i]])); 00339 } 00340 00341 void ComputeBlockFactor (FlatArray<int> block, int bw, FlatBandCholeskyFactors<TM> & inv) const; 00342 00344 virtual void MultAdd (TSCAL s, const BaseVector & x, BaseVector & y) const; 00345 00347 virtual void MultTransAdd (TSCAL s, const BaseVector & x, BaseVector & y) const; 00348 00350 virtual BaseVector * CreateVector () const 00351 { 00352 return mat.CreateVector(); 00353 } 00354 00355 00356 int Height() const { return mat.Height(); } 00357 virtual int VHeight() const { return mat.Height(); } 00358 int Width() const { return mat.Width(); } 00359 virtual int VWidth() const { return mat.Width(); } 00360 00361 00363 virtual void GSSmooth (BaseVector & x, const BaseVector & b, 00364 int steps = 1) const; 00365 00366 virtual void GSSmoothPartial (BaseVector & x, const BaseVector & b, 00367 BaseVector & y) const; 00368 00369 00370 00372 virtual void GSSmoothResiduum (BaseVector & x, const BaseVector & b, 00373 BaseVector & res, int steps = 1) const; 00374 00375 00377 virtual void GSSmoothBack (BaseVector & x, const BaseVector & b, 00378 int steps = 1) const; 00379 00380 00381 virtual void GSSmoothBackPartial (BaseVector & x, const BaseVector & b, 00382 BaseVector & y) const; 00383 00384 00385 00386 00387 00388 void SmoothBlock (int i, 00389 FlatVector<TVX> & x, 00390 // const FlatVector<TVX> & b, 00391 FlatVector<TVX> & y) const; 00392 00393 00395 virtual void GSSmoothNumbering (BaseVector & x, const BaseVector & b, 00396 const Array<int> & numbering, 00397 int forward = 1) const 00398 { 00399 ; 00400 } 00401 00402 /* 00403 virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const 00404 { 00405 // Array<BandCholeskyFactors<TM>*> invdiag; 00406 00407 00408 int nels = 0; 00409 for (int i = 0; i < blocktable.Size(); i++) 00410 if (invdiag[i]) 00411 { 00412 int s = invdiag[i]->Size(); 00413 int bw = invdiag[i]->BandWidth(); 00414 nels += s*bw - (bw * (bw-1)) / 2; 00415 nels += s; 00416 } 00417 mu.Append (new MemoryUsageStruct ("BlockJacSym Table", blocktable.NElements()*sizeof(int), 1)); 00418 mu.Append (new MemoryUsageStruct ("BlockJacSym", nels*sizeof(TM), 2*blocktable.Size())); 00419 } 00420 */ 00421 }; 00422 00423 00424 } 00425 00426 #endif