NGSolve  4.9
linalg/blockjacobi.hpp
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