SOFA API  99934cbc
Open source framework for multi-physics simuation
sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector > Class Template Reference

#include <BTDLinearSolver.h>

Inheritance diagram for sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >:

Detailed Description

template<class Matrix, class Vector>
class sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >

Linear system solver using Thomas Algorithm for Block Tridiagonal matrices

References: Conte, S.D., and deBoor, C. (1972). Elementary Numerical Analysis. McGraw-Hill, New York http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm http://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) http://www4.ncsu.edu/eos/users/w/white/www/white/ma580/chap2.5.PDF

Public Attributes

Data< boolf_verbose
 Dump system state at each iteration. More...
 
Data< boolproblem
 display debug informations about subpartSolve computation More...
 
Data< boolsubpartSolve
 Allows for the computation of a subpart of the system. More...
 
Data< boolverification
 verification of the subpartSolve More...
 
Data< booltest_perf
 verification of performance More...
 
helper::vector< SubMatrixalpha_inv
 
helper::vector< SubMatrixlambda
 
helper::vector< SubMatrixB
 
Matrix::InvMatrixType Minv
 
MysparseM H
 
MysparseMit H_it
 
Vector bwdContributionOnLH
 
Vector fwdContributionOnRH
 
Vector _rh_buf
 
SubVector _acc_rh_bloc
 
SubVector _acc_lh_bloc
 
Index current_bloc
 
Index first_block
 
std::vector< SubVectorVec_dRH
 
helper::vector< IndexnBlockComputedMinv
 
Vector Y
 
Data< intf_blockSize
 dimension of the blocks in the matrix More...
 

Public Member Functions

 SOFA_CLASS (SOFA_TEMPLATE2(BTDLinearSolver, Matrix, Vector), SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, Matrix, Vector))
 
void my_identity (SubMatrix &Id, const Index size_id)
 
void invert (SubMatrix &Inv, const BlocType &m)
 
void invert (Matrix &M) override
 
void computeMinvBlock (Index i, Index j)
 
double getMinvElement (Index i, Index j)
 
void solve (Matrix &, Vector &x, Vector &b) override
 Solve Mx=b. More...
 
bool addJMInvJt (defaulttype::BaseMatrix *result, defaulttype::BaseMatrix *J, double fact) override
 
void init_partial_solve () override
 Init the partial solve. More...
 
void partial_solve (ListIndex &Iout, ListIndex &Iin, bool NewIn) override
 *Matrix& M, Vector& result, Vector& rh, */ More...
 
void init_partial_inverse (const Index &nb, const Index &bsize)
 
template<class RMatrix , class JMatrix >
bool addJMInvJt (RMatrix &result, JMatrix &J, double fact)
 

Protected Member Functions

 BTDLinearSolver ()
 

Attribute details

template<class Matrix , class Vector >
SubVector sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::_acc_lh_bloc
template<class Matrix , class Vector >
SubVector sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::_acc_rh_bloc
template<class Matrix , class Vector >
Vector sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::_rh_buf
template<class Matrix , class Vector >
helper::vector<SubMatrix> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::alpha_inv
template<class Matrix , class Vector >
helper::vector<SubMatrix> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::B
template<class Matrix , class Vector >
Vector sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::bwdContributionOnLH
template<class Matrix , class Vector >
Index sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::current_bloc
template<class Matrix , class Vector >
Data<int> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::f_blockSize

dimension of the blocks in the matrix

template<class Matrix , class Vector >
Data<bool> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::f_verbose

Dump system state at each iteration.

template<class Matrix , class Vector >
Index sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::first_block
template<class Matrix , class Vector >
Vector sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::fwdContributionOnRH
template<class Matrix , class Vector >
MysparseM sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::H
template<class Matrix , class Vector >
MysparseMit sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::H_it
template<class Matrix , class Vector >
helper::vector<SubMatrix> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::lambda
template<class Matrix , class Vector >
Matrix::InvMatrixType sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::Minv
template<class Matrix , class Vector >
helper::vector<Index> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::nBlockComputedMinv
template<class Matrix , class Vector >
Data<bool> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::problem

display debug informations about subpartSolve computation

template<class Matrix , class Vector >
Data<bool> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::subpartSolve

Allows for the computation of a subpart of the system.

template<class Matrix , class Vector >
Data<bool> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::test_perf

verification of performance

template<class Matrix , class Vector >
std::vector<SubVector> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::Vec_dRH
template<class Matrix , class Vector >
Data<bool> sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::verification

verification of the subpartSolve

template<class Matrix , class Vector >
Vector sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::Y

Constructor details

template<class Matrix , class Vector >
sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::BTDLinearSolver ( )
inlineprotected

Function details

template<class Matrix , class Vector >
bool sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::addJMInvJt ( defaulttype::BaseMatrix result,
defaulttype::BaseMatrix J,
double  fact 
)
override

Multiply the inverse of the system matrix by the transpose of the given matrix, and multiply the result with the given matrix J

Parameters
resultthe variable where the result will be added
Jthe matrix J to use
Returns
false if the solver does not support this operation, of it the system matrix is not invertible
template<class Matrix , class Vector >
template<class RMatrix , class JMatrix >
bool sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::addJMInvJt ( RMatrix &  result,
JMatrix &  J,
double  fact 
)
template<class Matrix , class Vector >
void sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::computeMinvBlock ( Index  i,
Index  j 
)

[ inva0-l0(Minv10) Minv10t Minv20t Minv30t ] Minv = Uinv Linv = [ (Minv11)(-l0t) inva1-l1(Minv21) Minv21t Minv31t ] [ (Minv21)(-l0t) (Minv22)(-l1t) inva2-l2(Minv32) Minv32t ] [ (Minv31)(-l0t) (Minv32)(-l1t) (Minv33)(-l2t) inva3 ]

template<class Matrix , class Vector >
double sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::getMinvElement ( Index  i,
Index  j 
)
template<class Matrix , class Vector >
void sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::init_partial_inverse ( const Index nb,
const Index bsize 
)
template<class Matrix , class Vector >
void sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::init_partial_solve ( )
override

Init the partial solve.

template<class Matrix , class Vector >
void sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::invert ( SubMatrix Inv,
const BlocType m 
)
template<class Matrix , class Vector >
void sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::invert ( Matrix &  M)
override
template<class Matrix , class Vector >
void sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::my_identity ( SubMatrix Id,
const Index  size_id 
)

Factorize M

[ A0 C0 0  0  ]         [ a0 0  0  0  ] [ I  l0 0  0  ]

M = [ B1 A1 C1 0 ] = L U = [ B1 a1 0 0 ] [ 0 I l1 0 ] [ 0 B2 A2 C2 ] [ 0 B2 a2 0 ] [ 0 0 I l2 ] [ 0 0 B3 A3 ] [ 0 0 B3 a3 ] [ 0 0 0 I ] [ a0 a0l0 0 0 ] M = [ B1 B1l0+a1 a1l1 0 ] [ 0 B2 B2l1+a2 a2l2 ] [ 0 0 B3 B3l2+a3 ] L X = [ a0X0 B1X0+a1X1 B2X1+a2X2 B3X2+a3X3 ] [ inva0 0 0 0 ] Linv = [ -inva1B1inva0 inva1 0 0 ] [ inva2B2inva1B1inva0 -inva2B2inva1 inva2 0 ] [ -inva3B3inva2B2inva1B1inva0 inva3B3inva2B2inva1 -inva3B3inva2 inva3 ] U X = [ X0+l0X1 X1+l1X2 X2+l2X3 X3 ] Uinv = [ I -l0 l0l1 -l0l1l2 ] [ 0 I -l1 l1l2 ] [ 0 0 I -l2 ] [ 0 0 0 I ]

        [ (I+l0(I+l1(I+l2inva3B3)inva2B2)inva1B1)inva0 -l0(I+l1(I+l2inva3B3)inva2B2)inva1 l0l1(inva2+l2inva3B3inva2) -l0l1l2inva3 ]

Minv = Uinv Linv = [ -((I+l1(I+l2inva3B3)inva2B2)inva1B1)inva0 (I+l1(I+l2inva3B3)inva2B2)inva1 -l1(inva2+l2inva3B3inva2) l1l2inva3 ] [ (((I+l2inva3B3)inva2B2)inva1B1)inva0 -((I+l2inva3B3)inva2B2)inva1 inva2+l2inva3B3inva2 -l2inva3 ] [ -inva3B3inva2B2inva1B1inva0 inva3B3inva2B2inva1 -inva3B3inva2 inva3 ]

[ inva0-l0(Minv10) (-l0)(Minv11) (-l0)(Minv12) (-l0)(Minv13) ] Minv = Uinv Linv = [ (Minv11)(-B1inva0) inva1-l1(Minv21) (-l1)(Minv22) (-l1)(Minv23) ] [ (Minv21)(-B1inva0) (Minv22)(-B2inva1) inva2-l2(Minv32) (-l2)(Minv33) ] [ (Minv31)(-B1inva0) (Minv32)(-B2inva1) (Minv33)(-B3inva2) inva3 ]

if M is symmetric (Ai = Ait and Bi+1 = C1t) : li = invai*Ci = (invai)t*(Bi+1)t = (B(i+1)invai)t

               [ inva0-l0(Minv11)(-l0t)     Minv10t          Minv20t      Minv30t ]

Minv = Uinv Linv = [ (Minv11)(-l0t) inva1-l1(Minv22)(-l1t) Minv21t Minv31t ] [ (Minv21)(-l0t) (Minv22)(-l1t) inva2-l2(Minv33)(-l2t) Minv32t ] [ (Minv31)(-l0t) (Minv32)(-l1t) (Minv33)(-l2t) inva3 ]

template<class Matrix , class Vector >
void sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::partial_solve ( ListIndex Iout,
ListIndex Iin,
bool  NewIn 
)
override

*Matrix& M, Vector& result, Vector& rh, */

partial solve : b is accumulated db is a sparse vector that is added to b partial_x is a sparse vector (with sparse map given) that provide the result of M x = b+db Solve Mx=b

template<class Matrix , class Vector >
sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::SOFA_CLASS ( SOFA_TEMPLATE2(BTDLinearSolver< Matrix, Vector >, Matrix, Vector)  ,
SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, Matrix, Vector)   
)
template<class Matrix , class Vector >
void sofa::component::linearsolver::BTDLinearSolver< Matrix, Vector >::solve ( Matrix &  ,
Vector &  x,
Vector &  b 
)
override

Solve Mx=b.