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