1 /////////////////////////////////////////////////////////////////////////////// 2 // // 3 // The Template Matrix/Vector Library for C++ was created by Mike Jarvis // 4 // Copyright (C) 1998 - 2016 // 5 // All rights reserved // 6 // // 7 // The project is hosted at https://code.google.com/p/tmv-cpp/ // 8 // where you can find the current version and current documention. // 9 // // 10 // For concerns or problems with the software, Mike may be contacted at // 11 // mike_jarvis17 [at] gmail. // 12 // // 13 // This software is licensed under a FreeBSD license. The file // 14 // TMV_LICENSE should have bee included with this distribution. // 15 // It not, you can get a copy from https://code.google.com/p/tmv-cpp/. // 16 // // 17 // Essentially, you can use this software however you want provided that // 18 // you include the TMV_LICENSE file in any distribution that uses it. // 19 // // 20 /////////////////////////////////////////////////////////////////////////////// 21 22 23 //--------------------------------------------------------------------------- 24 // 25 // This file contains the code for doing division using 26 // Cholesky Decomposition. 27 // 28 // The algorithm is much like the LU decomposition, but we don't do 29 // any pivoting, and since the source matrix is symmetric, L = LT 30 // (or for Hermition, L = Lt). 31 // 32 33 34 #ifndef TMV_SymBandCHD_H 35 #define TMV_SymBandCHD_H 36 37 #include "tmv/TMV_SymDivider.h" 38 #include "tmv/TMV_BaseSymBandMatrix.h" 39 #include "tmv/TMV_BaseDiagMatrix.h" 40 41 namespace tmv { 42 43 // Decompose A into L*Lt. 44 // On output, L = A.lowerBand(). 45 template <typename T> 46 void CH_Decompose(SymBandMatrixView<T> A); 47 48 // Decompose Tridiagonal A into L D Lt 49 template <typename T> 50 void LDL_Decompose(SymBandMatrixView<T> A); 51 52 template <typename T, int A1> CH_Decompose(HermBandMatrix<T,A1> & A)53 inline void CH_Decompose(HermBandMatrix<T,A1>& A) 54 { CH_Decompose(A.view()); } 55 56 template <typename T, int A1> CH_Decompose(SymBandMatrix<T,A1> & A)57 inline void CH_Decompose(SymBandMatrix<T,A1>& A) 58 { CH_Decompose(A.view()); } 59 60 template <typename T, int A1> LDL_Decompose(HermBandMatrix<T,A1> & A)61 inline void LDL_Decompose(HermBandMatrix<T,A1>& A) 62 { LDL_Decompose(A.view()); } 63 64 template <typename T, int A1> LDL_Decompose(SymBandMatrix<T,A1> & A)65 inline void LDL_Decompose(SymBandMatrix<T,A1>& A) 66 { LDL_Decompose(A.view()); } 67 68 69 template <typename T> 70 class HermBandCHDiv : public SymDivider<T> 71 { 72 73 public : 74 75 // 76 // Constructors 77 // 78 79 HermBandCHDiv(const GenSymBandMatrix<T>& A, bool _inplace); 80 ~HermBandCHDiv(); 81 82 // 83 // Divider Versions of DivEq and Div 84 // 85 86 template <typename T1> 87 void doLDivEq(MatrixView<T1> m) const; 88 template <typename T1> 89 void doRDivEq(MatrixView<T1> m) const; 90 template <typename T1, typename T2> 91 void doLDiv(const GenMatrix<T1>& m1, MatrixView<T2> m0) const; 92 template <typename T1, typename T2> 93 void doRDiv(const GenMatrix<T1>& m1, MatrixView<T2> m0) const; 94 95 // 96 // Determinant, Inverse 97 // 98 99 T det() const; 100 TMV_RealType(T) logDet(T* sign) const; 101 template <typename T1> 102 void doMakeInverse(MatrixView<T1> minv) const; 103 template <typename T1> 104 void doMakeInverse(SymMatrixView<T1> minv) const; makeInverse(SymMatrixView<TMV_RealType (T)> sinv)105 inline void makeInverse(SymMatrixView<TMV_RealType(T)> sinv) const 106 { 107 TMVAssert(isReal(T())); 108 TMVAssert(sinv.size() == colsize()); 109 doMakeInverse(sinv); 110 } makeInverse(SymMatrixView<TMV_ComplexType (T)> sinv)111 inline void makeInverse(SymMatrixView<TMV_ComplexType(T)> sinv) const 112 { 113 TMVAssert(sinv.size() == colsize()); 114 TMVAssert(sinv.isherm()); 115 doMakeInverse(sinv); 116 } 117 void doMakeInverseATA(MatrixView<T> minv) const; 118 bool isSingular() const; 119 120 #include "tmv/TMV_AuxAllDiv.h" 121 122 // 123 // Access Decomposition 124 // 125 126 const BandMatrix<T> getL() const; 127 const DiagMatrix<T> getD() const; 128 const GenSymBandMatrix<T>& getLL() const; 129 130 bool checkDecomp(const BaseMatrix<T>& m, std::ostream* fout) const; 131 132 private : 133 struct HermBandCHDiv_Impl; 134 auto_ptr<HermBandCHDiv_Impl> pimpl; 135 136 ptrdiff_t colsize() const; 137 ptrdiff_t rowsize() const; 138 139 private : 140 141 HermBandCHDiv(const HermBandCHDiv<T>&); 142 HermBandCHDiv<T>& operator=(const HermBandCHDiv<T>&); 143 144 }; 145 146 } // namespace tmv 147 148 #endif 149