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