v0.14.0
Loading...
Searching...
No Matches
Functions
cholesky.hpp File Reference

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.

Functions

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...
 

Detailed Description

cholesky decomposition

-*- c++ -*-

Definition in file cholesky.hpp.

Function Documentation

◆ 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
MATRIXtype of input matrix
TRIAtype of lower triangular output matrix
Asquare symmetric positive definite input matrix (only the lower triangle is accessed)
Llower 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 //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}
static Index< 'L', 3 > L
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'k', 3 > k
constexpr AssemblyType A

◆ 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
MATRIXtype of matrix A
Ainput: square symmetric positive definite matrix (only the lower triangle is accessed)
Aoutput: 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 //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}

◆ 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
La triangular matrix
xinput: right hand side b; output: solution x
Examples
UnsaturatedFlow.hpp.

Definition at line 221 of file cholesky.hpp.

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}

◆ 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
MATRIXtype of matrix A
Ainput: square symmetric positive definite matrix (only the lower triangle is accessed)
Aoutput: 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 // 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}
FTensor::Index< 'i', SPACE_DIM > i
const double T