Namespaces
namespace	detail
namespace	details
Classes
class	BackwardGaussSeidelFunctor
class	ConjugateGradientFunctor
class	ForwardGaussSeidelFunctor
class	GMRESFunctor
class	IdentityPreconditioner
class	JacobiFunctor
class	SymmetricGaussSeidelFunctor
Functions
template<typename size_type , typename matrix_type , typename vector_type >
void	backsolve (size_type m, matrix_type const &A, vector_type &x, vector_type const &b)
template<typename T >
void	backward_gauss_seidel (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > &x, boost::numeric::ublas::vector< T > const &b)
template<typename T >
void	backward_gauss_seidel (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > &x, boost::numeric::ublas::vector< T > const &b, size_t const &max_iterations, size_t &iterations)
template<typename matrix_type , typename triangular_matrix_type >
size_t	cholesky_decompose (matrix_type const &A, triangular_matrix_type &L)
template<typename matrix_type >
size_t	cholesky_decompose (matrix_type &A)
template<typename matrix_type >
size_t	incomplete_cholesky_decompose (matrix_type &A)
template<typename triangular_matrix_type , typename vector_type >
void	cholesky_solve (triangular_matrix_type const &L, vector_type &x, ublas::lower)
template<typename T >
void	cholesky_solve (ublas::compressed_matrix< T > const &A, ublas::vector< T > &x, ublas::vector< T > const &b)
template<typename T >
void	cholesky_solve (ublas::matrix< T > const &A, ublas::vector< T > &x, ublas::vector< T > const &b)
template<typename matrix_type , typename vector_type >
void	conjugate_gradient (matrix_type const &A, vector_type &x, vector_type const &b, size_t const &max_iterations, typename vector_type::value_type const &epsilon, size_t &iterations)
template<typename matrix_type , typename vector_type >
void	conjugate_gradient (matrix_type const &A, vector_type &x, vector_type const &b, size_t &iterations)
template<typename matrix_type , typename vector_type >
void	conjugate_gradient (matrix_type const &A, vector_type &x, vector_type const &b)
template<typename value_type , typename matrix_type >
void	diag (ublas::vector< value_type > const &v, matrix_type &D)
template<typename T >
void	forward_gauss_seidel (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > &x, boost::numeric::ublas::vector< T > const &b)
template<typename T >
void	forward_gauss_seidel (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > &x, boost::numeric::ublas::vector< T > const &b, size_t const &max_iterations, size_t &iterations)
template<typename matrix_type >
void	generate_PD (size_t const &n, matrix_type &A)
template<typename matrix_type >
void	fast_generate_PD (size_t const &n, matrix_type &A)
template<typename matrix_type >
void	generate_PSD (size_t const &n, matrix_type &A, typename matrix_type::value_type const &fraction)
template<typename matrix_type >
void	generate_random (size_t const &m, size_t const &n, matrix_type &A)
template<typename value_type >
void	generate_random (size_t const &n, ublas::vector< value_type > &v)
template<typename matrix_type , typename vector_type , typename preconditioner_type >
void	gmres (matrix_type const &A, vector_type &x, vector_type const &b, typename vector_type::size_type const &max_iterations, typename vector_type::size_type const &max_restart_iterations, typename vector_type::value_type const &tolerance, typename vector_type::value_type &relative_residual_error, typename vector_type::size_type &used_inner_iterations, typename vector_type::size_type &used_outer_iterations, typename vector_type::size_type &status, preconditioner_type const &P)
template<typename matrix_type , typename vector_type >
void	gmres (matrix_type const &A, vector_type &x, vector_type const &b, typename vector_type::size_type const &max_iterations, typename vector_type::size_type const &max_restart_iterations, typename vector_type::value_type const &tolerance, typename vector_type::value_type &relative_residual_error, typename vector_type::size_type &used_inner_iterations, typename vector_type::size_type &used_outer_iterations, typename vector_type::size_type &status)
template<typename matrix_type , typename vector_type >
void	gmres (matrix_type const &A, vector_type &x, vector_type const &b, typename vector_type::size_type const &max_iterations, typename vector_type::size_type const &max_restart_iterations, typename vector_type::value_type const &tolerance)
template<typename matrix_type , typename vector_type >
void	gmres (matrix_type const &A, vector_type &x, vector_type const &b)
template<typename matrix_type >
void	gram_schmidt (matrix_type &A)
template<typename matrix_type >
bool	is_orthonormal (matrix_type const &A)
template<typename matrix_type >
bool	is_symmetric (matrix_type const &A)
template<typename T >
void	jacobi (ublas::compressed_matrix< T > const &A, ublas::vector< T > &x, ublas::vector< T > const &b)
template<typename T >
void	jacobi (ublas::compressed_matrix< T > const &A, ublas::vector< T > &x, ublas::vector< T > const &b, size_t const &max_iterations, T const &, size_t &iterations)
template<typename T >
void	jacobi (ublas::compressed_matrix< T > const &A, ublas::vector< T > &x, ublas::vector< T > const &b, size_t &iterations)
template<class matrix_type , class vector_type >
bool	lu (matrix_type const &A, vector_type &x, vector_type const &b)
template<class matrix_type >
bool	lu_invert (matrix_type const &A, matrix_type &invA)
template<class matrix_type , typename invert_functor >
void	moore_penrose_pseudoinverse (matrix_type const &A, matrix_type &invA, invert_functor const &invert)
template<class matrix_type >
void	lu_moore_penrose_pseudoinverse (matrix_type const &A, matrix_type &invA)
template<class matrix_type >
void	svd_moore_penrose_pseudoinverse (matrix_type const &A, matrix_type &invA)
template<typename T >
void	prod (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, boost::numeric::ublas::vector< T > &y)
template<typename T >
void	prod (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, T const &s, boost::numeric::ublas::vector< T > &y)
template<typename T >
void	prod_add (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, boost::numeric::ublas::vector< T > &y)
template<typename T >
void	prod_add (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, T const &s, boost::numeric::ublas::vector< T > &y)
template<typename T >
void	prod_add_rhs (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, boost::numeric::ublas::vector< T > const &b, boost::numeric::ublas::vector< T > &y)
template<typename T >
T	prod_row (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, typename boost::numeric::ublas::vector< T >::size_type i)
template<typename T >
void	prod_sub (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, boost::numeric::ublas::vector< T > &y)
template<typename T >
void	prod_sub_rhs (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, boost::numeric::ublas::vector< T > const &b, boost::numeric::ublas::vector< T > &y)
template<typename T >
void	prod_trans (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, boost::numeric::ublas::vector< T > &y)
template<typename T >
void	prod_trans (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, boost::numeric::ublas::vector< T > const &b, boost::numeric::ublas::vector< T > &y)
template<typename T >
void	residual (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > const &x, boost::numeric::ublas::vector< T > const &b, boost::numeric::ublas::vector< T > &r)
template<typename T , typename solver_function_type >
void	shur_system (ublas::compressed_matrix< T > &A_aa, ublas::compressed_matrix< T > const &A_ab, ublas::compressed_matrix< T > const &C, ublas::compressed_matrix< T > const &invD, ublas::vector< T > &rhs_a, ublas::vector< T > &rhs_b, ublas::vector< T > &dx_a, ublas::vector< T > &dx_b, solver_function_type const &solve)
template<typename ME >
void	svd (ublas::matrix_expression< ME > const &A, ublas::matrix< typename ME::value_type > &U, ublas::vector< typename ME::value_type > &s, ublas::matrix< typename ME::value_type > &V)
template<typename matrix_type , typename vector_type >
void	svd (matrix_type const &A, vector_type &x, vector_type const &b)
template<class matrix_type >
void	svd_invert (matrix_type const &A, matrix_type &invA)
template<typename T >
void	symmetric_gauss_seidel (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > &x, boost::numeric::ublas::vector< T > const &b)
template<typename T >
void	symmetric_gauss_seidel (boost::numeric::ublas::compressed_matrix< T > const &A, boost::numeric::ublas::vector< T > &x, boost::numeric::ublas::vector< T > const &b, size_t const &max_iterations, size_t &iterations)
template<class E , class T , class VE >
std::basic_ostream< E, T > &	operator<< (std::basic_ostream< E, T > &os, ublas::vector_expression< VE > const &v)
template<class T , class F >
std::ostream &	operator<< (std::ostream &os, ublas::compressed_matrix< T, F > const &A)
template<class E , class T , class ME >
std::basic_ostream< E, T > &	operator<< (std::basic_ostream< E, T > &os, ublas::matrix_expression< ME > const &m)
template<typename matrix_type >
bool	read_dlm_matrix (std::string const &filename, matrix_type &A)
template<typename vector_type >
bool	read_dlm_vector (std::string const &filename, vector_type &x)
template<typename matrix_type >
bool	write_dlm_matrix (std::string const &filename, matrix_type &A)
template<typename vector_type >
bool	write_dlm_vector (std::string const &filename, vector_type &x)

Function Documentation

template<typename size_type , typename matrix_type , typename vector_type >

void OpenTissue::math::big::backsolve	(	size_type	m,
		matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b
	)			`[inline]`

Solve Linear Upper Triangular Problem.

Parameters:

	m	The size of the problem.
	A	A square upper triangular matrix of dimension at least m m.
	x	Upon return this argument holds the result of x = A^{-1} b. Note x is of dimension b upon return.
	b	The right hand side vector of A x = b. Note b is of dimension at least m.

template<typename T >

void OpenTissue::math::big::backward_gauss_seidel	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > &	x,
		boost::numeric::ublas::vector< T > const &	b
	)			`[inline]`

Backward Gauss-Seidel Iteration. This function performs a singel backward Gauss-Seidel iterations, that is it given x^k it computes x^{k+1} = inv(D+U)(b - L x^k).

Parameters:

	A	The matrix.
	x	When called this parameter holds the current value of the iterate, upon return this arugment holds the new value of the iterate.
	b	The right hand side vector.

template<typename T >

void OpenTissue::math::big::backward_gauss_seidel	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > &	x,
		boost::numeric::ublas::vector< T > const &	b,
		size_t const &	max_iterations,
		size_t &	iterations
	)			`[inline]`

Backward Gauss-Seidel Iteration. This function performs a several backward Gauss-Seidel iterations, that is it given x^0 it computes x^{k+1} = inv(D+U)(b - L x^k) for solve k>0.

Parameters:

	A	The matrix.
	x	When called this parameter holds the current value of the iterate, upon return this arugment holds the new value of the iterate.
	b	The right hand side vector.
	max_iterations	The maximum number of iterations that can be taken.
	iterations	Upon return this argument holds the number of actual iterations that were used.

template<typename matrix_type , typename triangular_matrix_type >

size_t OpenTissue::math::big::cholesky_decompose	(	matrix_type const &	A,
		triangular_matrix_type &	L
	)			`[inline]`

Decompose the symmetric positive definit matrix A into product L L^T.

Parameters:

	A	Square symmetric positive definite input matrix (only the lower triangle is accessed).
	L	Upon return this argument holds the lower triangular output matrix.

Returns:: nonzero if decompositon fails (the value is 1 + the numer of the failing row)

template<typename matrix_type >

size_t OpenTissue::math::big::cholesky_decompose ( matrix_type & A ) [inline]

Decompose the symmetric positive definit matrix A into product L L^T.

Parameters:

A

Upon call this argument holds a square symmetric positive definite matrix (only the lower triangle is accessed). Upon return the lower triangle of A is replaced by the cholesky factor

Returns:: nonzero if decompositon fails (the value is 1 + the numer of the failing row)

template<typename T >

void OpenTissue::math::big::cholesky_solve	(	ublas::compressed_matrix< T > const &	A,
		ublas::vector< T > &	x,
		ublas::vector< T > const &	b
	)			`[inline]`

Solve Linear System A x = b using Cholesky Factorization.

Parameters:

	A	The matrix (must be postive definite),
	x	Upon return this argument contains the solution.
	b	The right hand side vector.

template<typename T >

void OpenTissue::math::big::cholesky_solve	(	ublas::matrix< T > const &	A,
		ublas::vector< T > &	x,
		ublas::vector< T > const &	b
	)			`[inline]`

Solve Linear System A x = b using Cholesky Factorization.

Parameters:

	A	The matrix (must be postive definite),
	x	Upon return this argument contains the solution.
	b	The right hand side vector.

template<typename triangular_matrix_type , typename vector_type >

void OpenTissue::math::big::cholesky_solve	(	triangular_matrix_type const &	L,
		vector_type &	x,
		ublas::lower
	)			`[inline]`

Solve system L L^T x = b inplace

Parameters:

	L	A triangular matrix.
	x	On invokation this argument holds the right hand side, upon return it holds the solution x.

template<typename matrix_type , typename vector_type >

void OpenTissue::math::big::conjugate_gradient	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b,
		size_t const &	max_iterations,
		typename vector_type::value_type const &	epsilon,
		size_t &	iterations
	)			`[inline]`

Conjugate Gradient Solver.

This implementation is based on code from Gunther Winkler. See: http://www-user.tu-chemnitz.de/~wgu/ublas/matrix_sparse_usage.html

Parameters:

	A	A symmetric positive definite matrix.
	x	Upon return this argument holds a the solution to the system A x = b
	b	The right hand side vector.
	max_iterations	The maximum number of iterates allowed.
	epsilon	The stopping threshold to be used.
	iterations	Upon return this argument holds the number of used iterations.

< If called multiple times, the iteration counter needs to be cleared!

< Residual vector.

< Search direction (gradient).

< The next Kyrlov subspace vector, A*x^{k+1}

template<typename matrix_type , typename vector_type >

void OpenTissue::math::big::conjugate_gradient	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b,
		size_t &	iterations
	)			`[inline]`

Conjugate Gradient Solver.

Parameters:

	A	A symmetric positive definite matrix.
	x	Upon return this argument holds a the solution to the system A x = b
	b	The right hand side vector.
	iterations	Upon return this argument holds the number of used iterations.

template<typename matrix_type , typename vector_type >

void OpenTissue::math::big::conjugate_gradient	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b
	)			`[inline]`

Conjugate Gradient Solver.

Parameters:

	A	A symmetric positive definite matrix.
	x	Upon return this argument holds a the solution to the system A x = b
	b	The right hand side vector.

template<typename value_type , typename matrix_type >

void OpenTissue::math::big::diag	(	ublas::vector< value_type > const &	v,
		matrix_type &	D
	)			`[inline]`

Generate Diagonal Matrix. This function is a convenience function that is usefull for quick-and-dirty initialization.

Parameters:

	v	A vector which holds the diagonal values of the resulting matrix.
	A	A square diagonal matrix holding the values of the vector along the diagonal.

template<typename matrix_type >

void OpenTissue::math::big::fast_generate_PD	(	size_t const &	n,
		matrix_type &	A
	)			`[inline]`

Generate Symmetric Positive Definite (PD) Matrix. This function is a convenience function that is usefull for quick-and-dirty initialization.

Parameters:

	n	The number of rows and columns in the resulting matrix.
	A	Upon return this argument holds an n-by-n PD matrix.

template<typename T >

void OpenTissue::math::big::forward_gauss_seidel	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > &	x,
		boost::numeric::ublas::vector< T > const &	b
	)			`[inline]`

Forward Gauss-Seidel Iteration. This function performs a single Gauss-Seidel iteration, that is it solves x^{k+1} = inv(D+L)(b - U x^k).

Parameters:

	A	The matrix.
	x	When called this parameter holds the current value of the iterate, upon return this arugment holds the new value of the iterate.
	b	The right hand side vector.

template<typename T >

void OpenTissue::math::big::forward_gauss_seidel	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > &	x,
		boost::numeric::ublas::vector< T > const &	b,
		size_t const &	max_iterations,
		size_t &	iterations
	)			`[inline]`

Forward Gauss-Seidel Iteration. This function performs a several Gauss-Seidel iterations, that is it given x^0 it computes x^{k+1} = inv(D+L)(b - U x^k) for solve k>0.

Parameters:

	A	The matrix.
	x	When called this parameter holds the current value of the iterate, upon return this arugment holds the new value of the iterate.
	b	The right hand side vector.
	max_iterations	The maximum number of iterations that can be taken.
	iterations	Upon return this argument holds the number of actual iterations that were used.

template<typename matrix_type >

void OpenTissue::math::big::generate_PD	(	size_t const &	n,
		matrix_type &	A
	)			`[inline]`

Generate Symmetric Positive Definite (PD) Matrix. This function is a convenience function that is usefull for quick-and-dirty initialization.

Parameters:

	n	The number of rows and columns in the resulting matrix.
	A	Upon return this argument holds an n-by-n PD matrix.

template<typename matrix_type >

void OpenTissue::math::big::generate_PSD	(	size_t const &	n,
		matrix_type &	A,
		typename matrix_type::value_type const &	fraction
	)			`[inline]`

Generate Symmetric Positive Semi-Definite (PSD) Matrix. This function is a convenience function that is usefull for quick-and-dirty initialization.

Notice, The eigenvalues of the resulting matrix is generated from random values between zero and one.

Parameters:

	n	The number of rows and columns in the resulting matrix.
	fraction	The (stocastic) fraction of zero-valued eigenvalues. If this is set to zero then the matrix is a positive definite matrix. If the fraction is set close to one then the number of zero-eigenvalues grow very large.
	A	Upon return this argument holds an n-by-n PSD matrix.

template<typename matrix_type >

void OpenTissue::math::big::generate_random	(	size_t const &	m,
		size_t const &	n,
		matrix_type &	A
	)			`[inline]`

Generate Random Matrix. This function is a convenience function that is usefull for quick-and-dirty initialization.

Parameters:

	m	The number of rows in the resulting matrix.
	n	The number of columns in the resulting matrix.
	A	Upon return this argument holds an m-by-n matrix with a random value between zero and one in each entry of the matrix.

template<typename value_type >

void OpenTissue::math::big::generate_random	(	size_t const &	n,
		ublas::vector< value_type > &	v
	)			`[inline]`

Generate Random Vector. This function is a convenience function that is usefull for quick-and-dirty initialization.

Parameters:

	n	The number of entries in the resulting vector.
	A	Upon return this argument holds an n-vector with a random value between zero and one in each entry.

template<typename matrix_type , typename vector_type , typename preconditioner_type >

void OpenTissue::math::big::gmres	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b,
		typename vector_type::size_type const &	max_iterations,
		typename vector_type::size_type const &	max_restart_iterations,
		typename vector_type::value_type const &	tolerance,
		typename vector_type::value_type &	relative_residual_error,
		typename vector_type::size_type &	used_inner_iterations,
		typename vector_type::size_type &	used_outer_iterations,
		typename vector_type::size_type &	status,
		preconditioner_type const &	P
	)			`[inline]`

GMRES Generalized Minimum Residual Method. GeneralizedMinimalResidualSolver solves the non-symmetric linear system Ax = b using the Generalized Minimal Residual method See: http://en.wikipedia.org/wiki/GMRES http://netlib2.cs.utk.edu/linalg/html_templates/report.html Saad's Iterative Methods for Sparse Linear Systems http://www-users.cs.umn.edu/~saad/books.html

We will try to use the same default values as the GMRES method implemented in Matlab. If a value is zero then it is assumed that the parameter were left unspecified by end-user.

From the Matlab help:

X = GMRES(A,B) attempts to solve the system of linear equations A*X = B for X. The N-by-N coefficient matrix A must be square and the right hand side column vector B must have length N. A may be a function returning A*X. This uses the unrestarted method with MIN(N,10) total iterations.

GMRES(A,B,RESTART) restarts the method every RESTART iterations. If RESTART is N or [] then GMRES uses the unrestarted method as above.

GMRES(A,B,RESTART,TOL) specifies the tolerance of the method. If TOL is [] then GMRES uses the default, 1e-6.

GMRES(A,B,RESTART,TOL,MAXIT) specifies the maximum number of outer iterations. Note: the total number of iterations is RESTART*MAXIT. If MAXIT is [] then GMRES uses the default, MIN(N/RESTART,10). If RESTART is N or [] then the total number of iterations is MAXIT.

Parameters:

	A	The system matrix.
	x	Upon return this argument holds the solution vector.
	b	The right hand side vector.
	max_itreations	The maximum number of outer iterations that can be used.
	max_restart_iterations	The restarting parameter for GMRES (see Saad's Iterative Methods for Sparse Linear Systems).
	tolerence	The threshold used in stop criteria
	relative_residual_error	Upon return this argument holds the accuracy of the last solution computed by invoking the solve methods.
	used_inner_iterations	Upon return this argument holds the number of used inner iterations.
	used_outer_iterations	Upon return this argument holds the number of used outer iterations.
	status	Upon return this argument holds the status of the GMRES method. status flag of GMRES, 0 : GMRES converged to the desired tolerance within maximum iterations. 1 : GMRES iterated maximum iterations but did not converge 2 : Preconditioner was ill-conditioned 3 : GMRES stagnated (two consecutive iterates were the same)
	P	A preconditioner.

template<typename matrix_type , typename vector_type >

void OpenTissue::math::big::gmres	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b,
		typename vector_type::size_type const &	max_iterations,
		typename vector_type::size_type const &	max_restart_iterations,
		typename vector_type::value_type const &	tolerance,
		typename vector_type::value_type &	relative_residual_error,
		typename vector_type::size_type &	used_inner_iterations,
		typename vector_type::size_type &	used_outer_iterations,
		typename vector_type::size_type &	status
	)			`[inline]`

template<typename matrix_type , typename vector_type >

void OpenTissue::math::big::gmres	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b,
		typename vector_type::size_type const &	max_iterations,
		typename vector_type::size_type const &	max_restart_iterations,
		typename vector_type::value_type const &	tolerance
	)			`[inline]`

template<typename matrix_type , typename vector_type >

void OpenTissue::math::big::gmres	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b
	)			`[inline]`

template<typename matrix_type >

void OpenTissue::math::big::gram_schmidt ( matrix_type & A ) [inline]

Gram-Schmidt Orthonormalization. This function orthonormalizes all the columns of the specified matrix.

Parameters:

A

Upon invokation this argument holds the matrix that should be orthonormalized. Upon return the argument holds the orthonormalized matrix. Thus the orthonormalization is done in-place.

template<typename matrix_type >

size_t OpenTissue::math::big::incomplete_cholesky_decompose ( matrix_type & A ) [inline]

Decompose the symmetric positive definit matrix A into product L L^T.

Parameters:

A

Upon invokation this argument holds a square symmetric positive definite matrix (only the lower triangle is accessed).Upon return the lower triangle of A is replaced by the cholesky factor.

Returns:: nonzero if decompositon fails (the value is 1 + the numer of the failing row)

template<typename matrix_type >

bool OpenTissue::math::big::is_orthonormal ( matrix_type const & A ) [inline]

Orthonormal Testing. This function is intended for debugging purposes it has not been implemented with focus on performance or accuracy. The function simply takes a brute force approach and tests if all the columns of the specified matrix are orthonormal.

Parameters:

A

The matrix that should be tested.

Returns:: If the specified matrix is orthonormal then the return value is true otherwise it is false.

template<typename matrix_type >

bool OpenTissue::math::big::is_symmetric ( matrix_type const & A ) [inline]

Symmetric Testing. This function is intended for debugging purposes it has not been implemented with focus on performance or accuracy. The function simply takes a brute force approach and tests if the specified matrix is symmetric.

Parameters:

A

The matrix that should be tested.

Returns:: If the specified matrix is symmetric then the return value is true otherwise it is false.

template<typename T >

void OpenTissue::math::big::jacobi	(	ublas::compressed_matrix< T > const &	A,
		ublas::vector< T > &	x,
		ublas::vector< T > const &	b
	)			`[inline]`

Jacobi Iteration Method. This function computes the result of a single iteration using the Jacobi method.

Parameters:

	A	A matrix, A = L + D + U, where D is diagonal, L strict lower triangular, and U strict upper triangular.
	x	Upon return this argument holds the value of x^{k+1} = D^{-1}(b - (L+U) x^k).
	b	The right hand side vector.

template<typename T >

void OpenTissue::math::big::jacobi	(	ublas::compressed_matrix< T > const &	A,
		ublas::vector< T > &	x,
		ublas::vector< T > const &	b,
		size_t const &	max_iterations,
		T const &	,
		size_t &	iterations
	)			`[inline]`

Jacobi Method. This function is capable of running several iterations of the Jacobi method.

Parameters:

	A	A matrix, A = L + D + U, where D is diagonal, L strict lower triangular, and U strict upper triangular.
	x	At invokation this argument holds the initial value of x^0, Upon return this argument holds the value of x^{k+1} = D^{-1}(b - (L+U) x^k) for some k>0.
	b	The right hand side vector.
	max_iterations	The maximum number of iterations that is allowed.
	epsilon	A stopping threshold (currently not used).
	iterations	Upon return this argument holds the number of used iterations.

template<typename T >

void OpenTissue::math::big::jacobi	(	ublas::compressed_matrix< T > const &	A,
		ublas::vector< T > &	x,
		ublas::vector< T > const &	b,
		size_t &	iterations
	)			`[inline]`

Jacobi Method. This function is capable of running several iterations of the Jacobi method. It uses a default setting of maximum 15 iterations and a stopping threshold of 10e-4.

This function is merely a convenience function with typical default settings.

Parameters:

	A	A matrix, A = L + D + U, where D is diagonal, L strict lower triangular, and U strict upper triangular.
	x	At invokation this argument holds the initial value of x^0, Upon return this argument holds the value of x^{k+1} = D^{-1}(b - (L+U) x^k) for some k>0.
	b	The right hand side vector.
	iterations	Upon return this argument holds the number of used iterations.

template<class matrix_type , class vector_type >

bool OpenTissue::math::big::lu	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b
	)			`[inline]`

Solve Linear System using LU decomposition.

Parameters:

	A	The matrix.
	x	Upon return this argument holds the solution
	b	The right hand side vector.

Returns:: If succesfull then the return value is true otherwise it is false.

template<class matrix_type >

bool OpenTissue::math::big::lu_invert	(	matrix_type const &	A,
		matrix_type &	invA
	)			`[inline]`

LU matrix inversion. This is a modified version of

http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?action=browse&diff=1&id=LU_Matrix_Inversion

Parameters:

	A	The matrix to be inverted.
	invA	Upon return this argument holds the inverted matrix.

Returns:: If succesfully inverted then true otherwise false.

template<class matrix_type >

void OpenTissue::math::big::lu_moore_penrose_pseudoinverse	(	matrix_type const &	A,
		matrix_type &	invA
	)			`[inline]`

Pseudoinverse using LU decompostion.

Parameters:

	A	Matrix to be inverted.
	invA	Upon return holds the pseudoinverse of A.

template<class matrix_type , typename invert_functor >

void OpenTissue::math::big::moore_penrose_pseudoinverse	(	matrix_type const &	A,
		matrix_type &	invA,
		invert_functor const &	invert
	)			`[inline]`

Moore-Penrose Pseudoinverse.

Parameters:

	A	The matrix to be inverted.
	invA	Upon return this argument holds the pseudoinverse matrix.
	invert	A functor object (or pointer) that specifies the matrix inversion method to be used.

< number of rows

< number of columns

template<class E , class T , class ME >

std::basic_ostream<E, T>& OpenTissue::math::big::operator<<	(	std::basic_ostream< E, T > &	os,
		ublas::matrix_expression< ME > const &	m
	)			`[inline]`

General Matrix Stream Output Operator.

Parameters:

	os	The stream that the matrix should be output to.
	m	The matrix that should be outputted.

Returns:: The stream that the matrix was output to.

template<class T , class F >

std::ostream& OpenTissue::math::big::operator<<	(	std::ostream &	os,
		ublas::compressed_matrix< T, F > const &	A
	)			`[inline]`

Compressed Matrix Stream Output Operator.

Parameters:

	os	The output stream the matrix should be output to.
	A	Compressed Matrix to output

Returns:: The stream that the matrix was written to.

template<class E , class T , class VE >

std::basic_ostream<E, T>& OpenTissue::math::big::operator<<	(	std::basic_ostream< E, T > &	os,
		ublas::vector_expression< VE > const &	v
	)			`[inline]`

Vector Stream Output Operator.

Parameters:

	os	The stream to output the vector to.
	v	A vector expression.

Returns:: The stream that the vector was written to.

template<typename T >

void OpenTissue::math::big::prod	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y = prod(A,x)

Parameters:

	A	A compressed matrix.
	x	A vector.
	y	Upon return this argument holds the result of A times x.

template<typename T >

void OpenTissue::math::big::prod	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		T const &	s,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y = prod(A,x)*s

Parameters:

	A	A compressed matrix.
	x	A vector.
	s	A scaling
	y	Upon return this argument holds the result of A times x times s.

template<typename T >

void OpenTissue::math::big::prod_add	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y += prod(A,x)

Parameters:

	A	A compressed matrix.
	x	A vector.
	y	Upon return this argument holds the result adding the value of A times x to the current value of the argument.

template<typename T >

void OpenTissue::math::big::prod_add	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		T const &	s,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y += prod(A,x)*s

Parameters:

	A	A compressed matrix.
	x	A vector.
	s	A scaling factor.
	y	Upon return this argument holds the result adding the value of A times x times s to the current value of the argument.

template<typename T >

void OpenTissue::math::big::prod_add_rhs	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		boost::numeric::ublas::vector< T > const &	b,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y = prod(A,x) + b

Parameters:

	A	The matrix.
	x	A vector to be multiplied with the matrix.
	b	The right hand side vector.
	y	Upon return this argument holds the value of the computation.

template<typename T >

T OpenTissue::math::big::prod_row	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		typename boost::numeric::ublas::vector< T >::size_type	i
	)			`[inline]`

Compute y_i = prod(row_i(A),x).

Parameters:

	A	A matrix
	x	A vector to be multipled with the matrix.
	i	The index of the i'th row value that should be computed.

Returns:: The value of the product of the i'th row with x-vector.

template<typename T >

void OpenTissue::math::big::prod_sub	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y -= prod(A,x)

Parameters:

	A	A compressed matrix.
	x	A vector.
	y	Upon return this argument holds the result subtracting the value of A times x to the current value of the argument.

template<typename T >

void OpenTissue::math::big::prod_sub_rhs	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		boost::numeric::ublas::vector< T > const &	b,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y = prod(A,x) - b

Parameters:

	A	The matrix.
	x	A vector to be multiplied with the matrix.
	b	The right hand side vector.
	y	Upon return this argument holds the value of the computation.

template<typename T >

void OpenTissue::math::big::prod_trans	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y = prod(trans(A),x)

Parameters:

	A	The matrix.
	x	A vector to be multiplied with the transpose of the specified matrix.
	y	Upon return this argument holds the result of the computations.

template<typename T >

void OpenTissue::math::big::prod_trans	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		boost::numeric::ublas::vector< T > const &	b,
		boost::numeric::ublas::vector< T > &	y
	)			`[inline]`

Compute y = prod(trans(A),x) + b

Parameters:

	A	The matrix.
	x	A vector to be multiplied with the transpose of the specified matrix.
	b	A vector to be added to the matrix vector product
	y	Upon return this argument holds the result of the computations.

template<typename matrix_type >

bool OpenTissue::math::big::read_dlm_matrix	(	std::string const &	filename,
		matrix_type &	A
	)			`[inline]`

Load a matrix from a ASCII text file (Fortran format) C++ Interface: ReadDLMmatrix

Description: Provides a way to read a matlab space delimited matrix

See Matlab Function reference, dlmread, dlmwrite.

Copyright: See COPYING file that comes with this distribution

Parameters:

	m	The matrix(vector) with compatible interface to read in.
	filename	The text filename, e.g., "A.dlm"

Returns:: Whether the read operation succeed.

template<typename vector_type >

bool OpenTissue::math::big::read_dlm_vector	(	std::string const &	filename,
		vector_type &	x
	)			`[inline]`

Load a vector from a ASCII text file (Fortran format) C++ Interface: ReadDLMmatrix

Description: Provides a way to read a matlab space delimited matrix

See Matlab Function reference, dlmread, dlmwrite.

Parameters:

	m	The vector with compatible interface to read in.
	filename	The text filename, e.g., "A.dlm"

Returns:: Whether the read operation succeed.

template<typename T >

void OpenTissue::math::big::residual	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > const &	x,
		boost::numeric::ublas::vector< T > const &	b,
		boost::numeric::ublas::vector< T > &	r
	)			`[inline]`

Compute r = b - prod(A,x).

Parameters:

	A	The matrix.
	x	The solution vector.
	b	The right hand side vector.
	r	Upon return this argument holds the residual value of the matrix equation A x = b.

template<typename T , typename solver_function_type >

void OpenTissue::math::big::shur_system	(	ublas::compressed_matrix< T > &	A_aa,
		ublas::compressed_matrix< T > const &	A_ab,
		ublas::compressed_matrix< T > const &	C,
		ublas::compressed_matrix< T > const &	invD,
		ublas::vector< T > &	rhs_a,
		ublas::vector< T > &	rhs_b,
		ublas::vector< T > &	dx_a,
		ublas::vector< T > &	dx_b,
		solver_function_type const &	solve
	)			`[inline]`

Compute Shur Equation. WARNING: This function performs an inplace solve, which means that some of the arguments are used as temporaries. That is the values of the arguments are overridden.

Given the partitioned matrix equation

| A_aa A_ab | �| x_a | = | rhs_a | | C D | | x_b | | rhs_b |

This method solves for the solution using a Shur Complement. From the bottom row we derive

x_b = inv(D)( rhs_b - C x_a ) (*1)

Substituting this into the top row and cleaning up one gets

A_aa x_a + A_ab inv(D)( rhs_b - C x_a ) = rhs_a A_aa x_a + A_ab inv(D) rhs_b - A_ab inv(D) C x_a = rhs_a (A_aa - A_ab inv(D) C) x_a = rhs_a - A_ab inv(D) rhs_b

Introducing the notation

M = (A_aa - A_ab inv(D) C) q = rhs_a - A_ab inv(D) rhs_b

We discover the Shur Matrix equation

M x_a = q (*2)

The idea is to solve for x_a using (*2) and then substituing into (*1) to find x_b.

Parameters:

	A_aa	The top-left most sub-block of the matrix equation. Overridden upon return.
	A_ab	The top-right most sub-block of the matrix equation.
	C	The bottom-left most sub-block of the matrix equation.
	invD	The inverse of the bottom-right most sub-block of the matrix equation.
	rhs_a	The top most sub-block of the right hand side vector. Overridden upon return.
	rhs_b	The bottom most sub-block of the right hand side vector. Overridden upon return.
	dx_a	Upon return this argument holds the top most sub-block of the solution vector.
	dx_b	Upon return this argument holds the bottom most sub-block of the solution vector.
	solve	A function that solves a linear system, A x = b. It must have the signature void ()(compressed_matrix const & A,vector & x, vector const & b).

template<typename matrix_type , typename vector_type >

void OpenTissue::math::big::svd	(	matrix_type const &	A,
		vector_type &	x,
		vector_type const &	b
	)			`[inline]`

Singular Value Decomposition Solver. First the function computes the singular value decomposition, A = U D VT, then it solves for the x-solution using this decomposition:

x = V inv(D) (UT (b))

Parameters:

	A	The matrix.
	x	Upon return this argument hold the solution to the linear system A x = b
	b	The right hand side vector.

template<typename ME >

void OpenTissue::math::big::svd	(	ublas::matrix_expression< ME > const &	A,
		ublas::matrix< typename ME::value_type > &	U,
		ublas::vector< typename ME::value_type > &	s,
		ublas::matrix< typename ME::value_type > &	V
	)			`[inline]`

Compute Singular Value Decomposition of a matrix.

A = U S VT

Parameters:

	A	The matrix.
	U	Upon return this matrix holds orthogonal columns.
	s	A vector containing the singular values.
	V	Upon return this matrix is an orthogonal matrix.

template<class matrix_type >

void OpenTissue::math::big::svd_invert	(	matrix_type const &	A,
		matrix_type &	invA
	)			`[inline]`

Matrix Inversion by Singular Value Decomposition.

Parameters:

	A	The matrix to be inverted.
	invA	Upon return this argument holds the inverted matrix.

template<class matrix_type >

void OpenTissue::math::big::svd_moore_penrose_pseudoinverse	(	matrix_type const &	A,
		matrix_type &	invA
	)			`[inline]`

Pseudoinverse using SVD decompostion.

Parameters:

	A	Matrix to be inverted.
	invA	Upon return holds the pseudoinverse of A.

template<typename T >

void OpenTissue::math::big::symmetric_gauss_seidel	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > &	x,
		boost::numeric::ublas::vector< T > const &	b,
		size_t const &	max_iterations,
		size_t &	iterations
	)			`[inline]`

Symmetric Gauss-Seidel Solver. This function is capable of performing several iterations of the symmetric Gauss Seidel iteration.

Parameters:

	A	A matrix.
	x	At invokation this argument holds the initial value of x^0, Upon return this argument holds the new value of iterate.
	b	The right hand side vector.
	max_iterations	The maximum number of iterations that is allowed.
	epsilon	A stopping threshold (currently not used).
	iterations	Upon return this argument holds the number of used iterations.

template<typename T >

void OpenTissue::math::big::symmetric_gauss_seidel	(	boost::numeric::ublas::compressed_matrix< T > const &	A,
		boost::numeric::ublas::vector< T > &	x,
		boost::numeric::ublas::vector< T > const &	b
	)			`[inline]`

Symmetric Gauss-Seidel Iteration. This function implements a single iteration of a symmetric Gauss Seidel method. That is it performs a forward Gauss-Seidel iteration followed by a backward Gauss-Seidel iteration.

Parameters:

	A	A matrix.
	x	Upon return this argument holds the new value of iterate x.
	b	The right hand side vector.

template<typename matrix_type >

bool OpenTissue::math::big::write_dlm_matrix	(	std::string const &	filename,
		matrix_type &	A
	)			`[inline]`

Save a matrix to a ASCII text file (Fortran format). See Matlab Function reference, dlmread, dlmwrite.

Parameters:

	A	The matrix that should be written.
	filename	The text filename, e.g., "A.dlm"

Returns:: Whether the read operation succeed.

template<typename vector_type >

bool OpenTissue::math::big::write_dlm_vector	(	std::string const &	filename,
		vector_type &	x
	)			`[inline]`

Write a vector to a ASCII text file (Fortran format). Provides a way to write a matlab space delimited matrix. See Matlab Function reference, dlmread, dlmwrite.

Parameters:

	filename	The text filename, e.g., "x.dlm"
	x	The vector that should be written.

Returns:: Whether the read operation succeed.

OpenTissue::math::big Namespace Reference

Namespaces

Classes

Functions

Function Documentation