Templated forward solver class. More...

#include <fwdsolver.h>

Inheritance diagram for TFwdSolver< T >:

Public Member Functions
	TFwdSolver (const QMMesh *mesh, LSOLVER linsolver, double tol=1e-10)
	Constructor. Creates a forward solver instance. More...

	TFwdSolver (const QMMesh mesh, const char solver, double tol=1e-10)
	Constructor. Creates a forward solver instance. More...

	TFwdSolver (const QMMesh *mesh, ParamParser &pp)
	Constructor. Creates a forward solver instance. More...

	~TFwdSolver ()
	Destructor. Destroys the forward solver instance.

void	SetDataScaling (DataScale scl)
	Set the default scaling for data returned by the projection methods. More...

DataScale	GetDataScaling () const
	Returns the current default setting for data scaling. More...

void	SetPrecon (PreconType type)
	Set the preconditioning method. More...

void	SetPhaseUnwrap (bool unwrap)

void	ReadParams (ParamParser &pp)
	Read solver parameters from a parameter file. More...

void	WriteParams (ParamParser &pp)
	Write the current solver settings to a parameter file. More...

void	SetLinSolver (const char *solver, double tol=1e-10)
	Set the solver type and tolerance. More...

LSOLVER	LinSolver () const
	Returns the current solver type. More...

double	GetLinSolverTol () const
	Returns the solver tolerance. More...

void	SetLinSolverMaxit (int maxit)
	Set max iteration count for iterative solver. More...

const QMMesh *	MeshPtr () const
	Returns a pointer to the FEM mesh associated with the forward solver instance. More...

void	Allocate ()
	Evaluates fill structure of system matrices and allocates dynamic memory for them. More...

void	AssembleSystemMatrix (const Solution &sol, double omega=0, bool elbasis=false)
	Construct the FEM system matrix from a set of parameter distributions for a given modulation frequency. More...

void	AssembleSystemMatrixComponent (const RVector &prm, int type)

void	AssembleMassMatrix (const Mesh *mesh=0)
	Construct the FEM mass matrix for time-dependent problems. More...

void	Reset (const Solution &sol, double omega=0, bool elbasis=false)
	Reset the forward solver by re-building the system matrices from a new set of parameters. More...

void	CalcField (const TVector< T > &qvec, TVector< T > &phi, IterativeSolverResult *res=0) const
	Calculate photon density field for a given source vector. More...

void	CalcFields (const TCompRowMatrix< T > &qvec, TVector< T > phi, IterativeSolverResult res=0) const
	Calculate photon density fields for all sources. More...

void	CalcFields (int q0, int q1, const TCompRowMatrix< T > &qvec, TVector< T > phi, IterativeSolverResult res=0) const
	Calculate photon density fields for a range of sources. More...

TVector< T >	ProjectSingle (int q, const TCompRowMatrix< T > &mvec, const TVector< T > &phi, DataScale scl=DATA_DEFAULT) const
	Return boundary data for a single source, given the corresponding photon density field. More...

TVector< T >	ProjectAll (const TCompRowMatrix< T > &mvec, const TVector< T > *dphi, DataScale scl=DATA_DEFAULT)
	Return boundary data for all sources and all detectors, given photon density fields for all sources. More...

TVector< T >	ProjectAll (const TCompRowMatrix< T > &qvec, const TCompRowMatrix< T > &mvec, const Solution &sol, double omega, DataScale scl=DATA_DEFAULT)
	Return boundary data for all sources and all detectors, given parameter distributions. More...

RVector	ProjectAll_real (const TCompRowMatrix< T > &mvec, const TVector< T > *dphi, DataScale scl=DATA_DEFAULT)
	Return boundary data for the complex case in a real vector. More...

FVector	ProjectAll_singlereal (const TCompRowMatrix< T > &mvec, const TVector< T > *dphi, DataScale scl=DATA_DEFAULT)

RVector	ProjectAll_real (const TCompRowMatrix< T > &qvec, const TCompRowMatrix< T > &mvec, const Solution &sol, double omega, DataScale scl=DATA_DEFAULT)
	Return boundary data for all sources and all detectors, given parameters defined via a Solution instance. More...

FVector	ProjectAll_singlereal (const TCompRowMatrix< T > &qvec, const TCompRowMatrix< T > &mvec, const Solution &sol, double omega, DataScale scl=DATA_DEFAULT)

Public Attributes
const QMMesh *	meshptr
	pointer to the associated FEM mesh

LSOLVER	solvertp
	linear solver type

IterativeMethod	method
	iterative solver method, if applicable

TCompRowMatrix< T > *	F
	FEM system matrix.

TCompRowMatrix< T > *	FL
	lower triangle of system matrix decomposition

TVector< T > *	Fd
	diagonal of Cholesky factorisation

TPreconditioner< T > *	precon
	preconditioner instance

TCompRowMatrix< T > *	B
	mass matrix; only used for time-domain problems

void *	SuperLU
	SuperLU solver engine.

double	iterative_tol
	iterative solver tolerance

int	iterative_maxit
	iterative solver max iterations (0 for auto)

PreconType	precontp
	preconditioner

Protected Member Functions
void	Setup ()

void	SetupType ()

void	DeleteType ()

RVector	UnfoldComplex (const TVector< T > &vec) const
	Unfold real and imaginary parts of a complex vector. More...

FVector	UnfoldSComplex (const TVector< T > &vec) const

Protected Attributes
DataScale	dscale
	default data scaling: DATA_LIN or DATA_LOG

TVector< T > *	pphi
	work buffer for field calculation

bool	unwrap_phase
	use phase unwrapping?

Detailed Description

template<class T>
class TFwdSolver< T >

Templated forward solver class.

TFwdSolver is a template class which encapsulates the FEM diffusion forward solver. It can be instantiated either as TFwdSolver<double> (or RFwdSolver for short) to describe real-valued problems (continuous-wave measurements), or as TFwdSolver<complex> (or CFwdSolver) for complex-valued problems (frequency domain problems).

The TFwdSolver class has an associated linear solver (defined in the constructor) for solving the linear FEM system. This can be either a direct solver (Cholesky for TFwdSolver<double> and LU for TFwdSolver<complex>), or an iterative method (CG, BiCG, BiCGSTAB, GMRES). For iterative methods, a tolerance limit is also required.

Each TFwdSolver instance has an FEM mesh associated with it. The mesh is defined by a call to Allocate, which also constructs the sparse system matrix structure based on the mesh connectivity graph.

The system matrix is built with a call to Reset for a given set of nodal optical coefficients. Subsequently, the photon density field (given a source vector) can be evaluated with CalcField and CalcFields, and the measurement values on the boundary can be obtained with Project and ProjectAll.

Constructor & Destructor Documentation

template<class T >

TFwdSolver< T >::TFwdSolver	(	const QMMesh *	mesh,
		LSOLVER	linsolver,
		double	tol = `1e-10`
	)

Constructor. Creates a forward solver instance.

Parameters

mesh	pointer to associated QMMesh object
linsolver	linear solver type (see Linear solver methods)
tol	linear solver tolerance (only used for iterative methods)

template<class T >

TFwdSolver< T >::TFwdSolver	(	const QMMesh *	mesh,
		const char *	solver,
		double	tol = `1e-10`
	)

Constructor. Creates a forward solver instance.

Parameters

mesh	pointer to associated QMMesh object
solver	linear solver type, provided as a string (see notes)
tol	linear solver tolerance (only used for iterative methods)

Note: For valid strings to define the linear solver, see SetLinSolver.

template<class T >

TFwdSolver< T >::TFwdSolver	(	const QMMesh *	mesh,
		ParamParser &	pp
	)

Constructor. Creates a forward solver instance.

Parameters

mesh	pointer to associated QMMesh object
pp	parser to read solver options from.

Note: For recognised items in the configuration file used by the parser, see the ReadParams method.

Member Function Documentation

template<class T >

void TFwdSolver< T >::Allocate ( )

Evaluates fill structure of system matrices and allocates dynamic memory for them.

Note: This function evaluates the fill structure of the sparse system matrix from the element connectivity graph of the associated mesh and allocates storage for it.; If a direct solver was selected, memory for the factorised matrix is also allocated, based on a symbolic Cholesky factorisation of the system matrix.; For iterative solvers, a simple diagonal precoditioner is generated.

template<class T >

void TFwdSolver< T >::AssembleMassMatrix ( const Mesh * mesh = 0 )

Construct the FEM mass matrix for time-dependent problems.

Parameters

mesh	mesh pointer

Note: If the provided mesh pointer is NULL, then the mesh from the preceeding call to Allocate is used instead. If no mesh has been assigned, the method fails.

See Also: Allocate, AssembleSystemMatrix

template<class T >

void TFwdSolver< T >::AssembleSystemMatrix	(	const Solution &	sol,
		double	omega = `0`,
		bool	elbasis = `false`
	)

Construct the FEM system matrix from a set of parameter distributions for a given modulation frequency.

Parameters

sol	Solution instance containing parameter vectors
omega	modulation frequency (cycles/ps)

Note: For template type double, parameter omega must be 0.; A call to Allocate must have preceeded the first call to AssembleSystemMatrix.

See Also: Reset

template<class T >

void TFwdSolver< T >::CalcField	(	const TVector< T > &	qvec,
		TVector< T > &	phi,
		IterativeSolverResult *	res = `0`
	)		const

Calculate photon density field for a given source vector.

Parameters

[in]	qvec	source vector (h-basis)
[out]	phi	photon density field (h-basis)
[out]	res	optional; if used, and if an iterative solver is used, the structure is filled with iteration count and relative error.

See Also: CalcFields

template<class T >

void TFwdSolver< T >::CalcFields	(	const TCompRowMatrix< T > &	qvec,
		TVector< T > *	phi,
		IterativeSolverResult *	res = `0`
	)		const

Calculate photon density fields for all sources.

Parameters

[in]	qvec	array of column source vectors
[out]	phi	array of photon density fields
[out]	res	optional; if used, and if an iterative solver is used, the structure is filled with the maximum iteration count and the maximum relative error for any source

Note: On call, phi must be an array of nq vectors, where nq is the number of rows in qvec.

See Also: CalcField

template<class T >

void TFwdSolver< T >::CalcFields	(	int	q0,
		int	q1,
		const TCompRowMatrix< T > &	qvec,
		TVector< T > *	phi,
		IterativeSolverResult *	res = `0`
	)		const

Calculate photon density fields for a range of sources.

Parameters

	q0	min source index
	q1	max source index + 1
	qvec	complete array of column source vectors
[out]	phi	array of photon density fields
[out]	res	optional; if used, and if an iterative solver is used, the structure is filled with the maximum iteration count and the maximum relative error for any source

Note: This function calculates the photon density fields of a sub-range q of sources (q0 <= q < q1); The returned array phi contains q1-q0 fields.

template<class T >

DataScale TFwdSolver< T >::GetDataScaling ( ) const

Returns the current default setting for data scaling.

Returns: current data scaling method: DATA_LIN (linear scaling) or DATA_LOG (logarithmic scaling).

See Also: SetDataScaling, ProjectAll, ProjectAll_real

template<class T >

double TFwdSolver< T >::GetLinSolverTol ( ) const

inline

Returns the solver tolerance.

Returns: Current solver tolerance value.

Note: The tolerance is only relevant for iterative solvers.

References TFwdSolver< T >::iterative_tol.

template<class T >

LSOLVER TFwdSolver< T >::LinSolver ( ) const

inline

Returns the current solver type.

Returns: Solver type enumeration value (see Linear solver methods).

References TFwdSolver< T >::solvertp.

template<class T >

const QMMesh* TFwdSolver< T >::MeshPtr ( ) const

inline

Returns a pointer to the FEM mesh associated with the forward solver instance.

Returns: Mesh pointer

References TFwdSolver< T >::meshptr.

template<class T >

TVector<T> TFwdSolver< T >::ProjectAll	(	const TCompRowMatrix< T > &	mvec,
		const TVector< T > *	dphi,
		DataScale	scl = `DATA_DEFAULT`
	)

Return boundary data for all sources and all detectors, given photon density fields for all sources.

Parameters

mvec	array of measurement column vectors
dphi	array of photon density fields for all sources (h-basis representation)
scl	output data scaling: linear or logarithmic

Returns: data vector (size nQM)

template<class T >

TVector<T> TFwdSolver< T >::ProjectAll	(	const TCompRowMatrix< T > &	qvec,
		const TCompRowMatrix< T > &	mvec,
		const Solution &	sol,
		double	omega,
		DataScale	scl = `DATA_DEFAULT`
	)

Return boundary data for all sources and all detectors, given parameter distributions.

Parameters

qvec	array of source column vectors
mvec	array of measurement column vectors
sol	solution containing optical parameters in h-basis
omega	modulation frequency (cycles/ps)
scl	output data scaling (linear or logarithmic)

template<class T >

RVector TFwdSolver< T >::ProjectAll_real	(	const TCompRowMatrix< T > &	mvec,
		const TVector< T > *	dphi,
		DataScale	scl = `DATA_DEFAULT`
	)

Return boundary data for the complex case in a real vector.

Parameters

mvec	array of measurement column vectors
dphi	array of photon density fields for all sources (h-basis representation)
scl	output data scaling: linear or logarithmic

Returns: data vector (size nQM*2)

Note: This function splits the complex data vector into real and imaginary parts, and returns them concatenated in a real vector, to be compliant with the standard TOAST data vector format.

template<class T >

RVector TFwdSolver< T >::ProjectAll_real	(	const TCompRowMatrix< T > &	qvec,
		const TCompRowMatrix< T > &	mvec,
		const Solution &	sol,
		double	omega,
		DataScale	scl = `DATA_DEFAULT`
	)

Return boundary data for all sources and all detectors, given parameters defined via a Solution instance.

Parameters

qvec	array of source column vectors
mvec	array of measurement column vectors
sol	solution instance (h-basis representation)
omega	modulation frequency (cycles/ps)
scl	output data scaling: linear or logarithmic

Returns: data vector (size nQM or nQM*2)

Note: This version always returns a real vector. For complex data, it returns the real and imaginary parts in sequential blocks of the vector.; This function constructs the system matrix from the solution, then calculates the fields for all sources, and obtains the projections by applying the measurement vectors.

template<class T >

TVector<T> TFwdSolver< T >::ProjectSingle	(	int	q,
		const TCompRowMatrix< T > &	mvec,
		const TVector< T > &	phi,
		DataScale	scl = `DATA_DEFAULT`
	)		const

Return boundary data for a single source, given the corresponding photon density field.

Parameters

q	source index (>= 0)
mvec	array of measurement column vectors
phi	photon density field for source q (h-basis representation)

Returns: data vector (size nQMref[q])

template<class T >

void TFwdSolver< T >::ReadParams ( ParamParser & pp )

Read solver parameters from a parameter file.

Parameters

pp	parser instance

Note

This method recognises the following items in a parameter file:

LINSOLVER	string	Linear solver method. The following choices are available: DIRECT (LU solver) CG, (conjugate gradient solver) BICG (biconjugate gradient solver), BICGSTAB (biconjugate gradient stabilised solver), GMRES (generalised minimum residual solver), GAUSSSEIDEL (Gauss-Seidel solver)
LINSOLVER_TOL	real	Tolerance value for iterative solvers

The LINSOLVER_TOL item is not required if LINSOLVER = DIRECT.

Any required items not found in the file are queried interactively from the user.

template<class T >

void TFwdSolver< T >::Reset	(	const Solution &	sol,
		double	omega = `0`,
		bool	elbasis = `false`
	)

Reset the forward solver by re-building the system matrices from a new set of parameters.

Parameters

sol	Solution instance containing the parameter vectors
omega	modulation frequency (cycles/ps)
elbasis	parameter vectors in sol are provided in piecewise constant element basis rather than nodal basis

Note: For template type double, parameter omega must be 0.; The system matrices must have been allocated (once) with a previous call to Allocate.; This function rebuilds the system matrix with a call to AssembleSystemMatrix. If a direct solver was selected, it also performs a Cholesky (double) or LU factorisation (complex).; If the mass matrix exists, it is rebuilt as well

Bug:

Resetting the mass matrix shouldn't be necessary, since it doesn't depend on the parameters.

See Also: AssembleSystemMatrix, AssembleMassMatrix

template<class T >

void TFwdSolver< T >::SetDataScaling ( DataScale scl )

Set the default scaling for data returned by the projection methods.

Parameters

scl	scaling method (DATA_LIN or DATA_LOG)
The	default scaling defined by this method is used by projection methods whenever they use a scaling parameter of DATA_DEFAULT.

Note: Before the first call to SetDataScaling, the default setting is DATA_LIN (linear data scaling)

See Also: GetDataScaling, ProjectAll, ProjectAll_real

template<class T >

void TFwdSolver< T >::SetLinSolver	(	const char *	solver,
		double	tol = `1e-10`
	)

Set the solver type and tolerance.

Parameters

solver	linear solver type, provided as a string (see notes)
tol	linear solver tolerance (only used for iterative methods)

Note

Valid options for the linear solver are:

DIRECT	Direct solver (LU)
CG	Conjugate gradient solver
BICG	Bi-conjugate gradient solver
BICGSTAB	Bi-conjugate gradient stabilised solver
GMRES	Generalised minimum residual solver

template<class T >

void TFwdSolver< T >::SetLinSolverMaxit ( int maxit )

inline

Set max iteration count for iterative solver.

Parameters

maxit max number of iterations (0 for auto)

Note: If the max iteration count is set to 0, it will be set to the dimension of the linear problem.

References TFwdSolver< T >::iterative_maxit.

template<class T >

void TFwdSolver< T >::SetPrecon ( PreconType type )

Set the preconditioning method.

Parameters

type	Preconditioner type.

Note: The default preconditioner is None.; If the preconditioner is constructed from a ParamParser object, the preconditioner may also be set via the LINSOLVER_PRECON tag in the parameter file.

See Also: ReadParams

template<class T >

RVector TFwdSolver< T >::UnfoldComplex ( const TVector< T > & vec ) const

protected

Unfold real and imaginary parts of a complex vector.

Parameters

vec	Input vector argument of template type

Returns: Unfolded real vector.

Note: For complex template types, this creates a real output vector of twice the size of the input vector, containing the real and imaginary parts of the vector.; For real template types, this simply returns the input vector.

template<class T >

void TFwdSolver< T >::WriteParams ( ParamParser & pp )

Write the current solver settings to a parameter file.

Parameters

pp	parser instance

Note: This method writes the LINSOLVER and (for iterative methods only) the LINSOLVER_TOL items to the parameter file. See ReadParams for details on the items.

The documentation for this class was generated from the following file:

/home/ucacmas/src/toast/toast_sf/toastpp-code/toast/trunk/src/libstoast/fwdsolver.h

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Detailed Description

template<class T> class TFwdSolver< T >

Constructor & Destructor Documentation

Member Function Documentation

template<class T>
class TFwdSolver< T >