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