Linear and nonlinear solvers

Linear \( Ax = b\) and non-linear \( f(x) = b\). More...

Collaboration diagram for Linear and nonlinear solvers:

Classes
struct	dg::AndersonAcceleration< ContainerType >
	Anderson Acceleration of Fixed Point/Richardson Iteration for the nonlinear equation \( f(x) = b\). More...
class	dg::BICGSTABl< ContainerType >
	Preconditioned BICGSTAB(l) method to solve \( Ax=b\). More...
class	dg::ChebyshevIteration< ContainerType >
	Preconditioned Chebyshev iteration for solving \( PAx=Pb\). More...
struct	dg::ChebyshevPreconditioner< Matrix, ContainerType >
	Chebyshev Polynomial Preconditioner \( C( A)\). More...
struct	dg::ModifiedChebyshevPreconditioner< Matrix, ContainerType >
	Approximate inverse Chebyshev Polynomial Preconditioner \( A^{-1} = \frac{c_0}{2} I + \sum_{k=1}^{r}c_kT_k( Z)\). More...
struct	dg::LeastSquaresPreconditioner< Matrix, InnerPreconditioner, ContainerType >
	Least Squares Polynomial Preconditioner \( M^{-1} s( AM^{-1})\). More...
class	dg::EVE< ContainerType >
	The Eigen-Value-Estimator (EVE) finds largest Eigenvalue of \( M^{-1}Ax = \lambda_\max x\). More...
struct	dg::DefaultSolver< ContainerType >
	PCG Solver class for solving \( (y-\alpha\hat I(t,y)) = \rho\). More...
class	dg::LGMRES< ContainerType >
	Functor class for the right preconditioned LGMRES method to solve \( Ax=b\). More...
class	dg::PCG< ContainerType >
	Preconditioned conjugate gradient method to solve \( Ax=b\). More...

Typedefs
template<class ContainerType >
using	dg::FixedPointIteration = AndersonAcceleration<ContainerType>
	If you are looking for fixed point iteration: it is a special case of Anderson Acceleration.

Functions
template<class T >
T	dg::create::lu_pivot (dg::SquareMatrix< T > &m, std::vector< unsigned > &p)
	LU Decomposition with partial pivoting.
template<class T >
dg::SquareMatrix< T >	dg::create::inverse (const dg::SquareMatrix< T > &in)
	Invert a square matrix.
template<class T >
void	dg::lu_solve (const dg::SquareMatrix< T > &lu, const std::vector< unsigned > &p, std::vector< T > &b)
	Solve the linear system with the LU decomposition.
template<class T >
dg::SquareMatrix< T >	dg::invert (const dg::SquareMatrix< T > &in)
	Compute inverse of square matrix (alias for dg::create::inverse)
template<typename UnaryOp >
int	dg::bisection1d (UnaryOp &&op, double &x_min, double &x_max, const double eps)
	Find a root of a 1d function \( f(x) = 0\).

Detailed Description

Linear \( Ax = b\) and non-linear \( f(x) = b\).

Typedef Documentation

FixedPointIteration

template<class ContainerType >

using dg::FixedPointIteration = AndersonAcceleration<ContainerType>

If you are looking for fixed point iteration: it is a special case of Anderson Acceleration.

Function Documentation

bisection1d()

template<typename UnaryOp >

int dg::bisection1d	(	UnaryOp &&	op,
		double &	x_min,
		double &	x_max,
		const double	eps )

Find a root of a 1d function \( f(x) = 0\).

in given boundaries using bisection

It is assumed that a sign change occurs at the root. Function jumps closer to the root by checking the sign.

double xmin = 0, xmax = 2;
int num = dg::bisection1d( [](double x){ return x*x - 2.;}, xmin, xmax, 1e-10);
// Num calls = 10* ln (10)/ln(2) = 33.2...
CHECK( num == 35);
CHECK( fabs( xmin - sqrt(2)) < 1e-9);
CHECK( fabs( xmax - sqrt(2)) < 1e-9);

Template Parameters

UnaryOp

unary function operator

Parameters

op	Function or Functor
x_min	left boundary, contains new left boundary on execution
x_max	right boundary, contains new right boundary on execution
eps	accuracy of the root finding. Algorithm successful if \( |x_{\max} - x_{\min}| < \varepsilon |x_{\max}| + \varepsilon \)

Returns

number of function calls to reach the desired accuracy

Exceptions

NoRoot1d	if no root lies between the given boundaries
std::runtime_error	if after 60 steps the accuracy wasn't reached

Note

If the root is found exactly then x_min = x_max

inverse()

template<class T >

dg::SquareMatrix< T > dg::create::inverse

(

const dg::SquareMatrix< T > &

in

)

Invert a square matrix.

using lu decomposition in combination with our accurate scalar products

For example

        // We can handle almost singular matrices:
        unsigned n = 10;
        dg::SquareMatrix<double> t( n, 1.);
        double d = 1.00 + 1e-3;
        for( unsigned u=0; u<n; u++)
            t(u,u) = d;
 
        // Determinant is ~ 1e-26
        auto t_inv = dg::invert( t);
 
        for( unsigned i=0; i<n; i++)
        for( unsigned j=0; j<n; j++)
        {
            double inv = i == j ? (d+n-2.0)/(d*(d+n-2.0) - (n-1.0))
                                : -1.0/(d*(d+n-2.0)-(n-1.0));
            INFO( "Item ("<<i<<" "<<j<<") Inv "<<t_inv(i,j)<<" Ana "
                          <<inv<<" diff "<<t_inv(i,j)-inv);
            CHECK( fabs((t_inv(i,j) - inv)/inv) < 1e-10);
        }

Template Parameters

T	value type

Parameters

in	input matrix

Returns

the inverse of in if it exists

Exceptions

std::runtime_error

if in is singular

Note

uses extended accuracy of dg::exblas which makes it quite robust against almost singular matrices

invert()

template<class T >

dg::SquareMatrix< T > dg::invert

(

const dg::SquareMatrix< T > &

in

)

Compute inverse of square matrix (alias for dg::create::inverse)

using lu decomposition in combination with our accurate scalar products

For example

        // We can handle almost singular matrices:
        unsigned n = 10;
        dg::SquareMatrix<double> t( n, 1.);
        double d = 1.00 + 1e-3;
        for( unsigned u=0; u<n; u++)
            t(u,u) = d;
 
        // Determinant is ~ 1e-26
        auto t_inv = dg::invert( t);
 
        for( unsigned i=0; i<n; i++)
        for( unsigned j=0; j<n; j++)
        {
            double inv = i == j ? (d+n-2.0)/(d*(d+n-2.0) - (n-1.0))
                                : -1.0/(d*(d+n-2.0)-(n-1.0));
            INFO( "Item ("<<i<<" "<<j<<") Inv "<<t_inv(i,j)<<" Ana "
                          <<inv<<" diff "<<t_inv(i,j)-inv);
            CHECK( fabs((t_inv(i,j) - inv)/inv) < 1e-10);
        }

Template Parameters

T	value type

Parameters

in	input matrix

Returns

the inverse of in if it exists

Exceptions

std::runtime_error

if in is singular

Note

uses extended accuracy of dg::exblas which makes it quite robust against almost singular matrices

lu_pivot()

template<class T >

T dg::create::lu_pivot	(	dg::SquareMatrix< T > &	m,
		std::vector< unsigned > &	p )

LU Decomposition with partial pivoting.

For example

        // Use lu_pivot to compute determinant
        unsigned n = 10;
        dg::SquareMatrix<double> t( n, 1.);
        double eps = 1e-3;
        for( unsigned u=0; u<n; u++)
            t(u,u) = (1. + eps);
 
        // Works for almost singular matrices
        std::vector<unsigned> p;
        double num = dg::create::lu_pivot( t, p);
        double det = pow( eps, n)*(1+n/eps); // ~ 1e-26
 
        INFO( "Det "<<num<<" Ana "<<det<<" diff "<<(num-det)/det);
        CHECK( fabs( num-det )/det < 1e-10);

Template Parameters

T	value type

Exceptions

std::runtime_error

if the matrix is singular

Parameters

m	contains lu decomposition of input on output (inplace transformation)
p	contains the pivot elements on output (will be resized)

Returns

determinant of m

Note

uses extended accuracy of dg::exblas which makes it quite robust against almost singular matrices

See also

dg::lu_solve

lu_solve()

template<class T >

void dg::lu_solve	(	const dg::SquareMatrix< T > &	lu,
		const std::vector< unsigned > &	p,
		std::vector< T > &	b )

Solve the linear system with the LU decomposition.

Template Parameters

T	value type

Parameters

lu	result of dg::create::lu_pivot
p	pivot vector from dg::create::lu_pivot
b	right hand side (contains solution on output)

See also

dg::create::lu_pivot