v0.13.1
cholesky.hpp
Go to the documentation of this file.
1/** -*- c++ -*- \file cholesky.hpp \brief cholesky decomposition */
2/*
3 - begin : 2005-08-24
4 - copyright : (C) 2005 by Gunter Winkler, Konstantin Kutzkow
5 - email : guwi17@gmx.de
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public
9 License as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11
12 This library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 Lesser General Public License for more details.
16
17 You should have received a copy of the GNU Lesser General Public
18 License along with this library; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20
21*/
22
23#ifndef _H_CHOLESKY_HPP_
24#define _H_CHOLESKY_HPP_
25
26
27#include <cassert>
28
29#include <boost/numeric/ublas/vector.hpp>
30#include <boost/numeric/ublas/vector_proxy.hpp>
31
32#include <boost/numeric/ublas/matrix.hpp>
33#include <boost/numeric/ublas/matrix_proxy.hpp>
34
35#include <boost/numeric/ublas/vector_expression.hpp>
36#include <boost/numeric/ublas/matrix_expression.hpp>
37
38#include <boost/numeric/ublas/triangular.hpp>
39
40//namespace ublas = boost::numeric::ublas;
41
42
43/** \brief decompose the symmetric positive definit matrix A into product L L^T.
44 *
45 * \param MATRIX type of input matrix
46 * \param TRIA type of lower triangular output matrix
47 * \param A square symmetric positive definite input matrix (only the lower triangle is accessed)
48 * \param L lower triangular output matrix
49 * \return nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
50 */
51template < class MATRIX, class TRIA >
52size_t cholesky_decompose(const MATRIX& A, TRIA& L)
53{
54 using namespace ublas;
55
56 //typedef typename MATRIX::value_type T;
57
58 assert( A.size1() == A.size2() );
59 assert( A.size1() == L.size1() );
60 assert( A.size2() == L.size2() );
61
62 const size_t n = A.size1();
63
64 for (size_t k=0 ; k < n; k++) {
65
66 double qL_kk = A(k,k) - inner_prod( project( row(L, k), range(0, k) ),
67 project( row(L, k), range(0, k) ) );
68
69 if (qL_kk <= 0) {
70 return 1 + k;
71 } else {
72 double L_kk = sqrt( qL_kk );
73 L(k,k) = L_kk;
74
75 matrix_column<TRIA> cLk(L, k);
76 project( cLk, range(k+1, n) )
77 = ( project( column(A, k), range(k+1, n) )
78 - prod( project(L, range(k+1, n), range(0, k)),
79 project(row(L, k), range(0, k) ) ) ) / L_kk;
80 }
81 }
82 return 0;
83}
84
85
86/** \brief decompose the symmetric positive definit matrix A into product L L^T.
87 *
88 * \param MATRIX type of matrix A
89 * \param A input: square symmetric positive definite matrix (only the lower triangle is accessed)
90 * \param A output: the lower triangle of A is replaced by the cholesky factor
91 * \return nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
92 */
93template < class MATRIX >
94size_t cholesky_decompose(MATRIX& A)
95{
96 using namespace ublas;
97
98 //typedef typename MATRIX::value_type T;
99
100 const MATRIX& A_c(A);
101
102 const size_t n = A.size1();
103
104 for (size_t k=0 ; k < n; k++) {
105
106 double qL_kk = A_c(k,k) - inner_prod( project( row(A_c, k), range(0, k) ),
107 project( row(A_c, k), range(0, k) ) );
108
109 if (qL_kk <= 0) {
110 return 1 + k;
111 } else {
112 double L_kk = sqrt( qL_kk );
113
114 matrix_column<MATRIX> cLk(A, k);
115 project( cLk, range(k+1, n) )
116 = ( project( column(A_c, k), range(k+1, n) )
117 - prod( project(A_c, range(k+1, n), range(0, k)),
118 project(row(A_c, k), range(0, k) ) ) ) / L_kk;
119 A(k,k) = L_kk;
120 }
121 }
122 return 0;
123}
124
125#if 0
126 using namespace ublas;
127
128 // Operations:
129 // n * (n - 1) / 2 + n = n * (n + 1) / 2 multiplications,
130 // n * (n - 1) / 2 additions
131
132 // Dense (proxy) case
133 template<class E1, class E2>
134 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
135 lower_tag, column_major_tag) {
136 std::cout << " is_lc ";
137 typedef typename E2::size_type size_type;
138 typedef typename E2::difference_type difference_type;
139 typedef typename E2::value_type value_type;
140
141 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
142 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
143 size_type size = e2 ().size ();
144 for (size_type n = 0; n < size; ++ n) {
145#ifndef BOOST_UBLAS_SINGULAR_CHECK
146 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
147#else
148 if (e1 () (n, n) == value_type/*zero*/())
149 singular ().raise ();
150#endif
151 value_type t = e2 () (n) / e1 () (n, n);
152 e2 () (n) = t;
153 if (t != value_type/*zero*/()) {
154 project( e2 (), range(n+1, size) )
155 .minus_assign( t * project( column( e1 (), n), range(n+1, size) ) );
156 }
157 }
158 }
159#endif
160
161
162
163
164
165/** \brief decompose the symmetric positive definit matrix A into product L L^T.
166 *
167 * \param MATRIX type of matrix A
168 * \param A input: square symmetric positive definite matrix (only the lower triangle is accessed)
169 * \param A output: the lower triangle of A is replaced by the cholesky factor
170 * \return nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
171 */
172template < class MATRIX >
174{
175 using namespace ublas;
176
177 typedef typename MATRIX::value_type T;
178
179 // read access to a const matrix is faster
180 const MATRIX& A_c(A);
181
182 const size_t n = A.size1();
183
184 for (size_t k=0 ; k < n; k++) {
185
186 double qL_kk = A_c(k,k) - inner_prod( project( row( A_c, k ), range(0, k) ),
187 project( row( A_c, k ), range(0, k) ) );
188
189 if (qL_kk <= 0) {
190 return 1 + k;
191 } else {
192 double L_kk = sqrt( qL_kk );
193
194 // aktualisieren
195 for (size_t i = k+1; i < A.size1(); ++i) {
196 T* Aik = A.find_element(i, k);
197
198 if (Aik != 0) {
199 *Aik = ( *Aik - inner_prod( project( row( A_c, k ), range(0, k) ),
200 project( row( A_c, i ), range(0, k) ) ) ) / L_kk;
201 }
202 }
203
204 A(k,k) = L_kk;
205 }
206 }
207
208 return 0;
209}
210
211
212
213
214/** \brief solve system L L^T x = b inplace
215 *
216 * \param L a triangular matrix
217 * \param x input: right hand side b; output: solution x
218 */
219template < class TRIA, class VEC >
220void
221cholesky_solve(const TRIA& L, VEC& x, ublas::lower)
222{
223 using namespace ublas;
224// ::inplace_solve(L, x, lower_tag(), typename TRIA::orientation_category () );
225 inplace_solve(L, x, lower_tag() );
226 inplace_solve(trans(L), x, upper_tag());
227}
228
229
230#endif
static Index< 'L', 3 > L
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t incomplete_cholesky_decompose(MATRIX &A)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:173
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'k', 3 > k
const double T
double A
constexpr double t
plate stiffness
Definition: plate.cpp:59