v0.14.0
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 }

◆ 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 }
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
convert.n
n
Definition: convert.py:82
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20