cholesky decomposition
More...
#include <cassert>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector_expression.hpp>
#include <boost/numeric/ublas/matrix_expression.hpp>
#include <boost/numeric/ublas/triangular.hpp>
Go to the source code of this file.
|
template<class MATRIX , class TRIA > |
size_t | cholesky_decompose (const MATRIX &A, TRIA &L) |
| decompose the symmetric positive definit matrix A into product L L^T. More...
|
|
template<class MATRIX > |
size_t | cholesky_decompose (MATRIX &A) |
| decompose the symmetric positive definit matrix A into product L L^T. More...
|
|
template<class MATRIX > |
size_t | incomplete_cholesky_decompose (MATRIX &A) |
| decompose the symmetric positive definit matrix A into product L L^T. More...
|
|
template<class TRIA , class VEC > |
void | cholesky_solve (const TRIA &L, VEC &x, ublas::lower) |
| solve system L L^T x = b inplace More...
|
|
cholesky decomposition
-*- c++ -*-
Definition in file cholesky.hpp.
◆ cholesky_decompose() [1/2]
template<class MATRIX , class TRIA >
size_t cholesky_decompose |
( |
const MATRIX & |
A, |
|
|
TRIA & |
L |
|
) |
| |
decompose the symmetric positive definit matrix A into product L L^T.
- Parameters
-
MATRIX | type of input matrix |
TRIA | type of lower triangular output matrix |
A | square symmetric positive definite input matrix (only the lower triangle is accessed) |
L | lower triangular output matrix |
- Returns
- nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
- Examples
- UnsaturatedFlow.hpp.
Definition at line 52 of file cholesky.hpp.
53{
54 using namespace ublas;
55
56
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) {
71 } else {
72 double L_kk = sqrt( qL_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}
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'k', 3 > k
◆ cholesky_decompose() [2/2]
template<class MATRIX >
size_t cholesky_decompose |
( |
MATRIX & |
A | ) |
|
decompose the symmetric positive definit matrix A into product L L^T.
- Parameters
-
MATRIX | type of matrix A |
A | input: square symmetric positive definite matrix (only the lower triangle is accessed) |
A | output: the lower triangle of A is replaced by the cholesky factor |
- Returns
- nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
Definition at line 94 of file cholesky.hpp.
95{
96 using namespace ublas;
97
98
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) {
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;
120 }
121 }
122 return 0;
123}
◆ cholesky_solve()
template<class TRIA , class VEC >
void cholesky_solve |
( |
const TRIA & |
L, |
|
|
VEC & |
x, |
|
|
ublas::lower |
|
|
) |
| |
solve system L L^T x = b inplace
- Parameters
-
L | a triangular matrix |
x | input: right hand side b; output: solution x |
- Examples
- UnsaturatedFlow.hpp.
Definition at line 221 of file cholesky.hpp.
222{
223 using namespace ublas;
224
225 inplace_solve(
L, x, lower_tag() );
226 inplace_solve(trans(
L), x, upper_tag());
227}
◆ incomplete_cholesky_decompose()
template<class MATRIX >
size_t incomplete_cholesky_decompose |
( |
MATRIX & |
A | ) |
|
decompose the symmetric positive definit matrix A into product L L^T.
- Parameters
-
MATRIX | type of matrix A |
A | input: square symmetric positive definite matrix (only the lower triangle is accessed) |
A | output: the lower triangle of A is replaced by the cholesky factor |
- Returns
- nonzero if decompositon fails (the value ist 1 + the numer of the failing row)
Definition at line 173 of file cholesky.hpp.
174{
175 using namespace ublas;
176
177 typedef typename MATRIX::value_type
T;
178
179
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) {
191 } else {
192 double L_kk = sqrt( qL_kk );
193
194
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
205 }
206 }
207
208 return 0;
209}
FTensor::Index< 'i', SPACE_DIM > i