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