1 | /*! |
---|
2 | * \file |
---|
3 | * \brief Matrices in decomposed forms (LDL', LU, UDU', etc). |
---|
4 | * \author Vaclav Smidl. |
---|
5 | * |
---|
6 | * ----------------------------------- |
---|
7 | * BDM++ - C++ library for Bayesian Decision Making under Uncertainty |
---|
8 | * |
---|
9 | * Using IT++ for numerical operations |
---|
10 | * ----------------------------------- |
---|
11 | */ |
---|
12 | |
---|
13 | #ifndef DC_H |
---|
14 | #define DC_H |
---|
15 | |
---|
16 | #include <itpp/itbase.h> |
---|
17 | |
---|
18 | using namespace itpp; |
---|
19 | |
---|
20 | /*! \brief Virtual class for representation of double symmetric matrices in square-root form. |
---|
21 | |
---|
22 | All operations defined on this class should be optimized for the chosed decomposition. |
---|
23 | */ |
---|
24 | class sqmat { |
---|
25 | public: |
---|
26 | /*! |
---|
27 | * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$. |
---|
28 | * @param v Vector forming the outer product to be added |
---|
29 | * @param w weight of updating; can be negative |
---|
30 | |
---|
31 | BLAS-2b operation. |
---|
32 | */ |
---|
33 | virtual void opupdt( const vec &v, double w ) =0; |
---|
34 | |
---|
35 | /*! \brief Conversion to full matrix. |
---|
36 | */ |
---|
37 | |
---|
38 | virtual mat to_mat() =0; |
---|
39 | |
---|
40 | /*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$ |
---|
41 | @param C multiplying matrix, |
---|
42 | @param trans if true, product $V = C'*V*C$ will be computed instead; |
---|
43 | */ |
---|
44 | virtual void mult_sym( const mat &C, bool trans=true ) =0; |
---|
45 | |
---|
46 | |
---|
47 | /*! |
---|
48 | \brief Logarithm of a determinant. |
---|
49 | |
---|
50 | */ |
---|
51 | virtual double logdet() =0; |
---|
52 | |
---|
53 | /*! |
---|
54 | \brief Multiplies square root of $V$ by vector $x$. |
---|
55 | |
---|
56 | Used e.g. in generating normal samples. |
---|
57 | */ |
---|
58 | virtual vec sqrt_mult(vec &v) =0; |
---|
59 | |
---|
60 | /*! |
---|
61 | \brief Evaluates quadratic form $x= v'*V*v$; |
---|
62 | |
---|
63 | */ |
---|
64 | virtual double qform(vec &v) =0; |
---|
65 | |
---|
66 | // //! easy version of the |
---|
67 | // sqmat inv(); |
---|
68 | |
---|
69 | //! Clearing matrix so that it corresponds to zeros. |
---|
70 | virtual void clear() =0; |
---|
71 | |
---|
72 | //! Reimplementing common functions of mat: cols(). |
---|
73 | virtual int cols() =0; |
---|
74 | |
---|
75 | //! Reimplementing common functions of mat: cols(). |
---|
76 | virtual int rows() =0; |
---|
77 | |
---|
78 | protected: |
---|
79 | int dim; |
---|
80 | }; |
---|
81 | |
---|
82 | |
---|
83 | /*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix. |
---|
84 | |
---|
85 | This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices; |
---|
86 | */ |
---|
87 | class fsqmat: sqmat { |
---|
88 | void opupdt( const vec &v, double w ); |
---|
89 | mat to_mat(); |
---|
90 | void mult_sym( const mat &C, bool trans=false ); |
---|
91 | void mult_sym( const mat &C, fsqmat &U, bool trans=false ); |
---|
92 | void inv(fsqmat &Inv); |
---|
93 | void clear(); |
---|
94 | |
---|
95 | |
---|
96 | //! Constructor |
---|
97 | fsqmat(const mat &M); |
---|
98 | |
---|
99 | /*! \brief Matrix inversion preserving the chosen form. |
---|
100 | |
---|
101 | @param Inv a space where the inverse is stored. |
---|
102 | |
---|
103 | */ |
---|
104 | virtual void inv(fsqmat* Inv); |
---|
105 | |
---|
106 | |
---|
107 | }; |
---|
108 | |
---|
109 | class ldmat: sqmat { |
---|
110 | public: |
---|
111 | |
---|
112 | //! Construct by copy of L and D. |
---|
113 | ldmat( const mat &L, const vec &D ); |
---|
114 | //! Construct by decomposition of full matrix V. |
---|
115 | ldmat( mat V ); |
---|
116 | ldmat (); |
---|
117 | |
---|
118 | // Reimplementation of compulsory operatios |
---|
119 | |
---|
120 | void opupdt( const vec &v, double w ); |
---|
121 | mat to_mat(); |
---|
122 | void mult_sym( const mat &C, bool trans=false ); |
---|
123 | void add ( const ldmat &ld2, double w=1.0 ); |
---|
124 | double logdet(); |
---|
125 | double qform(vec &v); |
---|
126 | // sqmat& operator -= ( const sqmat & ld2 ); |
---|
127 | void clear(); |
---|
128 | int cols(); |
---|
129 | int rows(); |
---|
130 | vec sqrt_mult(vec &v); |
---|
131 | |
---|
132 | /*! \brief Matrix inversion preserving the chosen form. |
---|
133 | |
---|
134 | @param Inv a space where the inverse is stored. |
---|
135 | |
---|
136 | */ |
---|
137 | virtual void inv(ldmat &Inv); |
---|
138 | |
---|
139 | /*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class. |
---|
140 | |
---|
141 | @param Inv a space where the inverse is stored. |
---|
142 | |
---|
143 | */ |
---|
144 | void mult_sym( const mat &C, ldmat &U, bool trans=false ); |
---|
145 | |
---|
146 | /*! \brief Transforms general $A'D0A$ into pure $L'DL$ |
---|
147 | |
---|
148 | The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$ |
---|
149 | |
---|
150 | @param A general matrix |
---|
151 | @param D0 general vector |
---|
152 | |
---|
153 | */ |
---|
154 | void ldform( mat &A, vec &D0 ); |
---|
155 | |
---|
156 | ldmat& operator += (const ldmat &ldA); |
---|
157 | ldmat& operator -= (const ldmat &ldA); |
---|
158 | ldmat& operator *= (double x); |
---|
159 | |
---|
160 | friend std::ostream &operator<< ( std::ostream &os, ldmat &sq ); |
---|
161 | |
---|
162 | protected: |
---|
163 | vec D; |
---|
164 | mat L; |
---|
165 | |
---|
166 | }; |
---|
167 | |
---|
168 | //////// Operations: |
---|
169 | |
---|
170 | inline ldmat& ldmat::operator += (const ldmat &ldA) {this->add(ldA);return *this;} |
---|
171 | inline ldmat& ldmat::operator -= (const ldmat &ldA) {this->add(ldA,-1.0);return *this;} |
---|
172 | inline int ldmat::cols(){return L.cols();} |
---|
173 | inline int ldmat::rows(){return L.rows();} |
---|
174 | |
---|
175 | #endif // DC_H |
---|