sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector > Class Template Reference

#include <BTDLinearSolver.h>

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

Detailed Description

template<class Matrix, class Vector>
class sofa::component::linearsolver::direct::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< bool >	d_verbose
	Dump system state at each iteration. More...
Data< bool >	d_problem
	display debug informations about subpartSolve computation More...
Data< bool >	d_subpartSolve
	Allows for the computation of a subpart of the system. More...
Data< bool >	d_verification
	verification of the subpartSolve More...
DeprecatedAndRemoved	d_blockSize
type::vector< SubMatrix >	alpha_inv
type::vector< SubMatrix >	lambda
type::vector< SubMatrix >	B
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< SubVector >	Vec_dRH
type::vector< Index >	nBlockComputedMinv
Vector	Y

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 (linearalgebra::BaseMatrix *result, linearalgebra::BaseMatrix *J, SReal 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)
void	parse (sofa::core::objectmodel::BaseObjectDescription *arg) override
Public Member Functions inherited from sofa::component::linearsolver::MatrixLinearSolver< Matrix, Vector, ThreadManager >
void	resetSystem ()
void	resizeSystem (Size n)
void	setSystemMBKMatrix (const core::MechanicalParams *mparams)
void	rebuildSystem (SReal, SReal)
void	setSystemLHVector (core::MultiVecDerivId v)
void	applySystemSolution ()
void	applyConstraintForce (const sofa::core::ConstraintParams *, sofa::core::MultiVecDerivId, const linearalgebra::BaseVector *)
void	computeResidual (const core::ExecParams *, linearalgebra::BaseVector *)
GraphScatteredVector *	createPersistentVector ()
linearalgebra::BaseMatrix *	getSystemBaseMatrix ()
linearalgebra::BaseVector *	getSystemRHBaseVector ()
linearalgebra::BaseVector *	getSystemLHBaseVector ()
void	checkLinearSystem ()
bool	addJMInvJtLocal (GraphScatteredMatrix *M, MatrixLinearSolver< GraphScatteredMatrix, GraphScatteredVector, NoThreadManager >::ResMatrixType *result, const MatrixLinearSolver< GraphScatteredMatrix, GraphScatteredVector, NoThreadManager >::JMatrixType *J, const SReal fact)
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	resetSystem ()
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	resizeSystem (Size)
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	setSystemMBKMatrix (const core::MechanicalParams *mparams)
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	rebuildSystem (SReal massFactor, SReal forceFactor)
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	setSystemLHVector (core::MultiVecDerivId v)
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	applySystemSolution ()
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API GraphScatteredVector *	createPersistentVector ()
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API linearalgebra::BaseMatrix *	getSystemBaseMatrix ()
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API linearalgebra::BaseVector *	getSystemRHBaseVector ()
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API linearalgebra::BaseVector *	getSystemLHBaseVector ()
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	applyConstraintForce (const sofa::core::ConstraintParams *, sofa::core::MultiVecDerivId, const linearalgebra::BaseVector *)
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	computeResidual (const core::ExecParams *params, linearalgebra::BaseVector *f)
SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API void	checkLinearSystem ()

Protected Member Functions
	BTDLinearSolver ()

Attribute details

_acc_lh_bloc

template<class Matrix , class Vector >

SubVector sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::_acc_lh_bloc

_acc_rh_bloc

template<class Matrix , class Vector >

SubVector sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::_acc_rh_bloc

_rh_buf

template<class Matrix , class Vector >

Vector sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::_rh_buf

alpha_inv

template<class Matrix , class Vector >

type::vector<SubMatrix> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::alpha_inv

B

template<class Matrix , class Vector >

type::vector<SubMatrix> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::B

bwdContributionOnLH

template<class Matrix , class Vector >

Vector sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::bwdContributionOnLH

current_bloc

template<class Matrix , class Vector >

Index sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::current_bloc

d_blockSize

template<class Matrix , class Vector >

DeprecatedAndRemoved sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::d_blockSize

d_problem

template<class Matrix , class Vector >

Data<bool> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::d_problem

display debug informations about subpartSolve computation

d_subpartSolve

template<class Matrix , class Vector >

Data<bool> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::d_subpartSolve

Allows for the computation of a subpart of the system.

d_verbose

template<class Matrix , class Vector >

Data<bool> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::d_verbose

Dump system state at each iteration.

d_verification

template<class Matrix , class Vector >

Data<bool> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::d_verification

verification of the subpartSolve

first_block

template<class Matrix , class Vector >

Index sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::first_block

fwdContributionOnRH

template<class Matrix , class Vector >

Vector sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::fwdContributionOnRH

H

template<class Matrix , class Vector >

MysparseM sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::H

H_it

template<class Matrix , class Vector >

MysparseMit sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::H_it

lambda

template<class Matrix , class Vector >

type::vector<SubMatrix> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::lambda

Minv

template<class Matrix , class Vector >

Matrix::InvMatrixType sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::Minv

nBlockComputedMinv

template<class Matrix , class Vector >

type::vector<Index> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::nBlockComputedMinv

Vec_dRH

template<class Matrix , class Vector >

std::vector<SubVector> sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::Vec_dRH

Y

template<class Matrix , class Vector >

Vector sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::Y

Constructor details

BTDLinearSolver()

template<class Matrix , class Vector >

sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::BTDLinearSolver ( )

inlineprotected

Function details

addJMInvJt() [1/2]

template<class Matrix , class Vector >

bool sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::addJMInvJt ( linearalgebra::BaseMatrix *  result, linearalgebra::BaseMatrix *  J, SReal  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

result	the variable where the result will be added
J	the matrix J to use

Returns

false if the solver does not support this operation, of it the system matrix is not invertible

addJMInvJt() [2/2]

template<class Matrix , class Vector >

template<class RMatrix , class JMatrix >

bool sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::addJMInvJt	(	RMatrix &	result,
		JMatrix &	J,
		double	fact
	)

computeMinvBlock()

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::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 ]

getMinvElement()

template<class Matrix , class Vector >

double sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::getMinvElement	(	Index	i,
		Index	j
	)

init_partial_inverse()

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::init_partial_inverse	(	const Index &	nb,
		const Index &	bsize
	)

init_partial_solve()

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::init_partial_solve

override

Init the partial solve.

invert() [1/2]

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::invert ( Matrix &  M )

override

invert() [2/2]

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::invert	(	SubMatrix &	Inv,
		const BlocType &	m
	)

my_identity()

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::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 ]

parse()

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::parse ( sofa::core::objectmodel::BaseObjectDescription *  arg )

inlineoverride

partial_solve()

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::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

SOFA_CLASS()

template<class Matrix , class Vector >

sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::SOFA_CLASS	(	SOFA_TEMPLATE2(BTDLinearSolver< Matrix, Vector >, Matrix, Vector)	,
		SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, Matrix, Vector)
	)

solve()

template<class Matrix , class Vector >

void sofa::component::linearsolver::direct::BTDLinearSolver< Matrix, Vector >::solve ( Matrix &  , Vector &  x, Vector &  b  )

override

Solve Mx=b.