NGSolve  4.9
basiclinalg/cholesky.hpp
00001 #ifndef FILE_CHOLESKY
00002 #define FILE_CHOLESKY
00003 
00004 /****************************************************************************/
00005 /* File:   cholesky.hpp                                                     */
00006 /* Author: Joachim Schoeberl                                                */
00007 /* Date:   25. Mar. 2000, 16. June 2002                                     */
00008 /****************************************************************************/
00009 
00010 namespace ngbla
00011 {
00012 
00017   template <class T>
00018   class FlatCholeskyFactors
00019   {
00020   protected:
00022     int n;
00024     T * lfact;
00026     T * diag;
00027   public:
00028     typedef typename mat_traits<T>::TV_COL TV;
00030     FlatCholeskyFactors (const FlatMatrix<T> & a, T * data)
00031     {
00032       diag = data;
00033       Factor (a);
00034     }
00035 
00037     FlatCholeskyFactors (const FlatMatrix<T> & a, LocalHeap & lh)
00038     {
00039       diag = (T*)lh.Alloc(sizeof(T)*RequiredMem(a.Height()));
00040       Factor (a);
00041     }
00042 
00044     NGS_DLL_HEADER void Factor (const FlatMatrix<T> & a);
00046     NGS_DLL_HEADER void Mult (const FlatVector<TV> & x, FlatVector<TV> & y) const;
00048     NGS_DLL_HEADER ostream & Print (ostream & ost) const;
00049 
00050 
00052     static int RequiredMem (int n)
00053     { return n*(n+1)/2; }
00054 
00055   private:
00057     T * PRow (int i) const { return lfact + (i*(i-1)) / 2; }
00058   };
00059 
00060 
00062   template<typename T>
00063   inline std::ostream & operator<< (std::ostream & s, const FlatCholeskyFactors<T> & m)
00064   {
00065     m.Print (s);
00066     return s;
00067   }
00068 
00069 
00070 
00071   template <class T>
00072   class CholeskyFactors : public FlatCholeskyFactors<T>
00073   {
00074   public:
00076     CholeskyFactors (const FlatMatrix<T> & a)
00077       : FlatCholeskyFactors<T> (a, new T[this->RequiredMem(a.Height())])
00078     { ; }
00080     ~CholeskyFactors ()
00081     {
00082       delete [] this->diag;
00083     }
00084   };
00085 
00086 }
00087 
00088 #endif