NGSolve
4.9
|
00001 #ifndef FILE_SUPERLUINVERSE 00002 #define FILE_SUPERLUINVERSE 00003 00004 /* *************************************************************************/ 00005 /* File: superluinverse.hpp */ 00006 /* Author: Florian Bachinger */ 00007 /* Date: Feb. 04 */ 00008 /* *************************************************************************/ 00009 00010 namespace ngla 00011 { 00012 00013 /* 00014 interface to the sparse direct solver SuperLU 00015 */ 00016 00017 00018 #ifdef USE_SUPERLU 00019 00020 namespace double_superlu 00021 { 00022 //#include "superlu/dsp_defs.h" 00023 //#include "superlu/util.h" 00024 extern "C" 00025 { 00026 00027 #include <slu_ddefs.h> 00028 #include <slu_util.h> 00029 #include <supermatrix.h> 00030 } 00031 00032 } 00033 namespace complex_superlu 00034 { 00035 extern "C" 00036 { 00037 using namespace double_superlu; 00038 #include <slu_zdefs.h> 00039 #include <slu_util.h> 00040 #include <supermatrix.h> 00041 #include <slu_scomplex.h> 00042 } 00043 00044 } 00045 00046 using namespace double_superlu; 00047 using namespace complex_superlu; 00048 00049 #endif 00050 00052 00053 00054 template<class TM, 00055 class TV_ROW = typename mat_traits<TM>::TV_ROW, 00056 class TV_COL = typename mat_traits<TM>::TV_COL> 00057 class SuperLUInverse : public BaseMatrix 00058 { 00059 #ifdef USE_SUPERLU 00060 SuperMatrix A, B, L, U; 00061 00062 SuperLUStat_t stat; 00063 #endif 00064 int height, nze, entrysize; 00065 00066 int * perm_c, * perm_r; 00067 int * colstart, * indices; 00068 typename mat_traits<TM>::TSCAL * matrix; 00069 00070 int symmetric, iscomplex; 00071 00072 const BitArray * inner; 00073 const Array<int> * cluster; 00074 00075 public: 00076 typedef typename mat_traits<TM>::TV_COL TV; 00077 typedef typename mat_traits<TM>::TV_ROW TVX; 00078 typedef typename mat_traits<TM>::TSCAL TSCAL; 00079 00081 SuperLUInverse (const SparseMatrix<TM,TV_ROW,TV_COL> & a, 00082 const BitArray * ainner = NULL, 00083 const Array<int> * acluster = NULL, 00084 int symmetric = 0); 00085 // /// 00086 // SuperLUInverse (const SparseMatrix<TM> & a, 00087 // const BitArray * ainner = NULL, 00088 // const Array<int> * acluster = NULL, 00089 // int symmetric = 0); 00090 00092 SuperLUInverse (const Array<int> & aorder, 00093 const Array<CliqueEl*> & cliques, 00094 const Array<MDOVertex> & vertices, 00095 int symmetric = 0); 00097 ~SuperLUInverse (); 00099 int VHeight() const { return height; } 00101 int VWidth() const { return height; } 00103 void Allocate (const Array<int> & aorder, 00104 const Array<CliqueEl*> & cliques, 00105 const Array<MDOVertex> & vertices); 00107 void Factor (const int * blocknr); 00109 void FactorNew (const SparseMatrix<TM> & a); 00111 virtual void Mult (const BaseVector & x, BaseVector & y) const; 00112 00114 virtual ostream & Print (ostream & ost) const; 00115 00116 virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const 00117 { 00118 mu.Append (new MemoryUsageStruct ("SuperLUInverse", nze*sizeof(TM), 1)); 00119 } 00120 00121 /* /// 00122 void Set (int i, int j, const TM & val); 00124 const TM & Get (int i, int j) const; 00125 */ 00126 virtual BaseVector * CreateVector () const 00127 { 00128 return new VVector<TV> (height); 00129 } 00130 }; 00131 00132 } 00133 00134 #endif