NGSolve
4.9
|
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