157 #include <boost/numeric/ublas/matrix_sparse.hpp>
159 #include <feel/feelcore/feel.hpp>
160 #include <feel/feelcore/traits.hpp>
166 namespace ublas= boost::numeric::ublas;
168 template<
typename T =
double>
172 typedef T value_type;
173 typedef ublas::compressed_matrix<value_type, ublas::column_major> matrix_type;
175 ldl( matrix_type& A )
178 Ap( A.index1_data().begin() ),
179 Ai( A.index2_data().begin() ),
180 Ax( A.value_data().begin() ),
197 void lsolve ( ublas::vector<value_type>& X );
199 void dsolve ( ublas::vector<value_type>& X );
201 void ltsolve ( ublas::vector<value_type>& X );
206 void solve ( ublas::vector<value_type>& X )
213 void perm ( ublas::vector<value_type>& X, ublas::vector<value_type>
const& B );
214 void permt ( ublas::vector<value_type>& X, ublas::vector<value_type>
const& B );
217 bool valid_matrix ();
227 ublas::vector<size_type> Lp;
229 ublas::vector<int64_type> Parent;
232 ublas::vector<size_type> Lnz;
233 ublas::vector<size_type> Li;
234 ublas::vector<value_type> Lx;
235 ublas::vector<value_type> D;
236 ublas::vector<value_type> Y;
239 ublas::vector<size_type> Flag;
240 ublas::vector<size_type> Pattern;
277 void ldl<T>::symbolic()
284 for ( k = 0 ; k < n ; k++ )
290 for ( k = 0 ; k < n ; k++ )
296 kk = ( P ) ? ( P [k] ) : ( k ) ;
299 for ( p = Ap [kk] ; p < p2 ; p++ )
302 i = ( Pinv ) ? ( Pinv [Ai [p]] ) : ( Ai [p] ) ;
307 for ( ; Flag [i] != k ; i = Parent [i] )
310 if ( Parent [i] == -1 ) Parent [i] = k ;
322 for ( k = 0 ; k < n ; k++ )
324 Lp [k+1] = Lp [k] + Lnz [k] ;
346 value_type yi, l_ki ;
350 for ( k = 0 ; k < n ; k++ )
357 kk = ( P ) ? ( P [k] ) : ( k ) ;
360 for ( p = Ap [kk] ; p < p2 ; p++ )
362 i = ( Pinv ) ? ( Pinv [Ai [p]] ) : ( Ai [p] ) ;
368 for ( len = 0 ; Flag [i] != k ; i = Parent [i] )
370 Pattern [len++] = i ;
374 while ( len > 0 ) Pattern [--top] = Pattern [--len] ;
382 for ( ; top < n ; top++ )
387 p2 = Lp [i] + Lnz [i] ;
389 for ( p = Lp [i] ; p < p2 ; p++ )
391 Y [Li [p]] -= Lx [p] * yi ;
401 if ( D [k] == 0.0 )
return ( k ) ;
412 void ldl<T>::lsolve( ublas::vector<value_type>& X )
416 for ( j = 0 ; j < n ; j++ )
420 for ( p = Lp [j] ; p < p2 ; p++ )
422 X [Li [p]] -= Lx [p] * X [j] ;
432 void ldl<T>::dsolve( ublas::vector<value_type>& X )
436 for ( j = 0 ; j < n ; j++ )
447 void ldl<T>::ltsolve( ublas::vector<value_type>& X )
451 for (
int j = n-1 ; j >= 0 ; j-- )
455 for ( p = Lp [j] ; p < p2 ; p++ )
457 X [j] -= Lx [p] * X [Li [p]] ;
467 void ldl<T>::perm( ublas::vector<value_type> & X,
468 ublas::vector<value_type>
const& B )
472 for ( j = 0 ; j < n ; j++ )
483 void ldl<T>::permt( ublas::vector<value_type> & X,
484 ublas::vector<value_type>
const& B )
488 for ( j = 0 ; j < n ; j++ )
501 bool ldl<T>::valid_perm()
505 if ( Flag.size() == 0 )
515 for ( j = 0 ; j < n ; j++ )
520 for ( k = 0 ; k < n ; k++ )
524 if ( j >= n || Flag [j] != 0 )
547 bool ldl<T>::valid_matrix()
556 for ( j = 0 ; j < n ; j++ )
558 if ( Ap [j] > Ap [j+1] )
564 for ( p = 0 ; p < Ap [n] ; p++ )
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319