MathLib::LinAlg Namespace Reference

Detailed Description

Some general linear algebra functionality.

By using the provided functions linear algebra capabilities can be used for different matrix and vector types in a way that is agnostic towards the specific type used.

For documentation, refer to that of the templated method. All specializations or overload must behave in the same way.

Functions
void	setLocalAccessibleVector (PETScVector const &x)
void	set (PETScVector &x, PetscScalar const a)
void	copy (PETScVector const &x, PETScVector &y)
void	scale (PETScVector &x, PetscScalar const a)
void	aypx (PETScVector &y, PetscScalar const a, PETScVector const &x)
void	axpy (PETScVector &y, PetscScalar const a, PETScVector const &x)
void	axpby (PETScVector &y, PetscScalar const a, PetscScalar const b, PETScVector const &x)
template<>
void	componentwiseDivide (PETScVector &w, PETScVector const &x, PETScVector const &y)
template<>
double	norm1 (PETScVector const &x)
template<>
double	norm2 (PETScVector const &x)
template<>
double	normMax (PETScVector const &x)
void	copy (PETScMatrix const &A, PETScMatrix &B)
void	scale (PETScMatrix &A, PetscScalar const a)
void	aypx (PETScMatrix &Y, PetscScalar const a, PETScMatrix const &X)
void	axpy (PETScMatrix &Y, PetscScalar const a, PETScMatrix const &X)
void	matMult (PETScMatrix const &A, PETScVector const &x, PETScVector &y)
void	matMultAdd (PETScMatrix const &A, PETScVector const &v1, PETScVector const &v2, PETScVector &v3)
void	linearSysNormalize (PETScMatrix const &, PETScMatrix &, PETScVector const &, PETScVector &)
void	finalizeAssembly (PETScMatrix &A)
void	finalizeAssembly (PETScVector &x)
template<typename MatrixOrVector >
void	copy (MatrixOrVector const &x, MatrixOrVector &y)
	Copies x to y.
template<typename MatrixOrVector >
void	scale (MatrixOrVector &x, double const a)
	Scales x by a.
template<typename MatrixOrVector >
void	aypx (MatrixOrVector &y, double const a, MatrixOrVector const &x)
	Computes $y = a \cdot y + x$.
template<typename MatrixOrVector >
void	axpy (MatrixOrVector &y, double const a, MatrixOrVector const &x)
	Computes $y = a \cdot x + y$.
template<typename MatrixOrVector >
void	axpby (MatrixOrVector &y, double const a, double const b, MatrixOrVector const &x)
	Computes $y = a \cdot x + b \cdot y$.
template<typename MatrixOrVector >
void	componentwiseDivide (MatrixOrVector &w, MatrixOrVector const &x, MatrixOrVector const &y)
	Computes $w = x/y$ componentwise.
template<typename MatrixOrVector >
double	norm1 (MatrixOrVector const &x)
	Computes the Manhattan norm of x.
template<typename MatrixOrVector >
double	norm2 (MatrixOrVector const &x)
	Computes the Euclidean norm of x.
template<typename MatrixOrVector >
double	normMax (MatrixOrVector const &x)
	Computes the maximum norm of x.
template<typename MatrixOrVector >
double	norm (MatrixOrVector const &x, MathLib::VecNormType type)
template<typename Matrix >
void	finalizeAssembly (Matrix &)
template<typename Matrix , typename Vector >
void	matMult (Matrix const &A, Vector const &x, Vector &y)
template<typename Matrix , typename Vector >
void	matMultAdd (Matrix const &A, Vector const &v1, Vector const &v2, Vector &v3)
template<typename VectorType >
double	computeRelativeNorm (VectorType const &x, VectorType const &y, MathLib::VecNormType norm_type)

Function Documentation

axpby() [1/2]

template<typename MatrixOrVector >

void MathLib::LinAlg::axpby	(	MatrixOrVector &	y,
		double const	a,
		double const	b,
		MatrixOrVector const &	x )

Computes $y = a \cdot x + b \cdot y$.

Definition at line 72 of file LinAlg.h.

{
    y = a * x + b * y;
}

axpby() [2/2]

void MathLib::LinAlg::axpby	(	PETScVector &	y,
		PetscScalar const	a,
		PetscScalar const	b,
		PETScVector const &	x )

Definition at line 64 of file LinAlg.cpp.

{
    // TODO check sizes
    VecAXPBY(y.getRawVector(), a, b, x.getRawVector());
}

References MathLib::PETScVector::getRawVector().

axpy() [1/3]

template<typename MatrixOrVector >

void MathLib::LinAlg::axpy	(	MatrixOrVector &	y,
		double const	a,
		MatrixOrVector const &	x )

Computes $y = a \cdot x + y$.

Definition at line 65 of file LinAlg.h.

{
    y += a * x;
}

axpy() [2/3]

void MathLib::LinAlg::axpy	(	PETScMatrix &	Y,
		PetscScalar const	a,
		PETScMatrix const &	X )

Definition at line 139 of file LinAlg.cpp.

{
    // TODO check sizes
    // TODO sparsity pattern, currently they are assumed to be different (slow)
    MatAXPY(Y.getRawMatrix(), a, X.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
}

References MathLib::PETScMatrix::getRawMatrix().

axpy() [3/3]

void MathLib::LinAlg::axpy	(	PETScVector &	y,
		PetscScalar const	a,
		PETScVector const &	x )

Definition at line 57 of file LinAlg.cpp.

{
    // TODO check sizes
    VecAXPY(y.getRawVector(), a, x.getRawVector());
}

References MathLib::PETScVector::getRawVector().

Referenced by ProcessLib::AssemblyMixin< Process >::assembleWithJacobian(), NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::calculateNonEquilibriumInitialResiduum(), NumLib::StaggeredCoupling::checkCouplingConvergence(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeA(), computeRelativeNorm(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeResidual(), NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::solve(), and NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::solve().

aypx() [1/3]

template<typename MatrixOrVector >

void MathLib::LinAlg::aypx	(	MatrixOrVector &	y,
		double const	a,
		MatrixOrVector const &	x )

Computes $y = a \cdot y + x$.

Definition at line 58 of file LinAlg.h.

{
    y = a * y + x;
}

aypx() [2/3]

void MathLib::LinAlg::aypx	(	PETScMatrix &	Y,
		PetscScalar const	a,
		PETScMatrix const &	X )

Definition at line 131 of file LinAlg.cpp.

{
    // TODO check sizes
    // TODO sparsity pattern, currently they are assumed to be different (slow)
    MatAYPX(Y.getRawMatrix(), a, X.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
}

References MathLib::PETScMatrix::getRawMatrix().

aypx() [3/3]

void MathLib::LinAlg::aypx	(	PETScVector &	y,
		PetscScalar const	a,
		PETScVector const &	x )

Definition at line 50 of file LinAlg.cpp.

{
    // TODO check sizes
    VecAYPX(y.getRawVector(), a, x.getRawVector());
}

References MathLib::PETScVector::getRawVector().

Referenced by ProcessLib::ComponentTransport::ComponentTransportProcess::solveReactionEquation().

componentwiseDivide() [1/2]

template<typename MatrixOrVector >

void MathLib::LinAlg::componentwiseDivide	(	MatrixOrVector &	w,
		MatrixOrVector const &	x,
		MatrixOrVector const &	y )

Computes $w = x/y$ componentwise.

componentwiseDivide() [2/2]

template<>

void MathLib::LinAlg::componentwiseDivide	(	PETScVector &	w,
		PETScVector const &	x,
		PETScVector const &	y )

Definition at line 81 of file LinAlg.cpp.

{
    VecPointwiseDivide(w.getRawVector(), x.getRawVector(), y.getRawVector());
}

References MathLib::PETScVector::getRawVector().

Referenced by NumLib::LocalLinearLeastSquaresExtrapolator::extrapolate().

computeRelativeNorm()

template<typename VectorType >

double MathLib::LinAlg::computeRelativeNorm	(	VectorType const &	x,
		VectorType const &	y,
		MathLib::VecNormType	norm_type )

Compute the relative norm between two vectors by $\|x-y\|/\|x\|$.

Parameters

x	Vector x
y	Vector y
norm_type	The norm type of global vector

Returns

$\|x-y\|/\|x\|$.

Definition at line 299 of file LinAlg.h.

{
    if (norm_type == MathLib::VecNormType::INVALID)
    {
        OGS_FATAL("An invalid norm type has been passed");
    }
 
    // Stores \f$diff = x-y\f$.
    VectorType diff;
    MathLib::LinAlg::copy(x, diff);
    MathLib::LinAlg::axpy(diff, -1.0, y);
 
    const double norm_diff = MathLib::LinAlg::norm(diff, norm_type);
 
    const double norm_x = MathLib::LinAlg::norm(x, norm_type);
    if (norm_x > std::numeric_limits<double>::epsilon())
    {
        return norm_diff / norm_x;
    }
 
    // Both of norm_x and norm_diff are close to zero
    if (norm_diff < std::numeric_limits<double>::epsilon())
    {
        return 1.0;
    }
 
    // Only norm_x is close to zero
    return norm_diff / std::numeric_limits<double>::epsilon();
}

References axpy(), copy(), MathLib::INVALID, norm(), and OGS_FATAL.

Referenced by ProcessLib::TimeLoop::computeTimeStepping().

copy() [1/3]

template<typename MatrixOrVector >

void MathLib::LinAlg::copy	(	MatrixOrVector const &	x,
		MatrixOrVector &	y )

Copies x to y.

Definition at line 44 of file LinAlg.h.

{
    y = x;
}

copy() [2/3]

void MathLib::LinAlg::copy	(	PETScMatrix const &	A,
		PETScMatrix &	B )

Definition at line 119 of file LinAlg.cpp.

{
    B = A;
}

copy() [3/3]

void MathLib::LinAlg::copy	(	PETScVector const &	x,
		PETScVector &	y )

Definition at line 37 of file LinAlg.cpp.

{
    if (!y.getRawVector())
        y.shallowCopy(x);
    VecCopy(x.getRawVector(), y.getRawVector());
}

References MathLib::PETScVector::getRawVector(), and MathLib::PETScVector::shallowCopy().

Referenced by ProcessLib::AssembledMatrixCache::assemble(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeA(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeJacobian(), computeRelativeNorm(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeResidual(), ProcessLib::TimeLoop::computeTimeStepping(), NumLib::TimeDiscretizedODESystem< ODESystemTag::FirstOrderImplicitQuasilinear, NonlinearSolverTag::Newton >::getResidual(), NumLib::SimpleMatrixVectorProvider::getVector(), NumLib::SimpleMatrixVectorProvider::getVector(), NumLib::BackwardEuler::getWeightedOldX(), NumLib::StaggeredCoupling::initializeCoupledSolutions(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::normalizeAandRhs(), ProcessLib::Process::setInitialConditions(), NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::solve(), NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::solve(), NumLib::PETScNonlinearSolver::solve(), ProcessLib::ComponentTransport::ComponentTransportProcess::solveReactionEquation(), ProcessLib::HMPhaseField::HMPhaseFieldProcess< DisplacementDim >::updateConstraints(), ProcessLib::PhaseField::PhaseFieldProcess< DisplacementDim >::updateConstraints(), and NumLib::StaggeredCoupling::updatePreviousSolution().

finalizeAssembly() [1/3]

template<typename Matrix >

void MathLib::LinAlg::finalizeAssembly

(

Matrix &

)

finalizeAssembly() [2/3]

void MathLib::LinAlg::finalizeAssembly

(

PETScMatrix &

A

)

Definition at line 198 of file LinAlg.cpp.

{
    A.finalizeAssembly(MAT_FINAL_ASSEMBLY);
}

References MathLib::PETScMatrix::finalizeAssembly().

Referenced by detail::applyKnownSolutions(), NumLib::TimeDiscretizedODESystem< ODESystemTag::FirstOrderImplicitQuasilinear, NonlinearSolverTag::Newton >::assemble(), NumLib::TimeDiscretizedODESystem< ODESystemTag::FirstOrderImplicitQuasilinear, NonlinearSolverTag::Picard >::assemble(), ProcessLib::AssemblyMixin< Process >::assembleWithJacobian(), NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::calculateNonEquilibriumInitialResiduum(), NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::calculateNonEquilibriumInitialResiduum(), NumLib::LocalLinearLeastSquaresExtrapolator::calculateResiduals(), NumLib::LocalLinearLeastSquaresExtrapolator::extrapolate(), ProcessLib::Process::setInitialConditions(), and ProcessLib::SmallDeformation::writeMaterialForces().

finalizeAssembly() [3/3]

void MathLib::LinAlg::finalizeAssembly

(

PETScVector &

x

)

Definition at line 203 of file LinAlg.cpp.

{
    x.finalizeAssembly();
}

References MathLib::PETScVector::finalizeAssembly().

linearSysNormalize()

void MathLib::LinAlg::linearSysNormalize	(	PETScMatrix const &	,
		PETScMatrix &	,
		PETScVector const &	,
		PETScVector &	)

Definition at line 170 of file LinAlg.cpp.

{
    // The following block is deactivated, because there is no tests yet for the
    // normalization operation in PETSc. This will be a task for later.
    /*
    assert(&A != &new_A);
    assert(&b != &new_b);
 
    PetscInt n_rows(0);
    PetscInt n_cols(0);
    MatGetSize(A.getRawMatrix(), &n_rows, &n_cols);
    // only when A matrix is not square
    if (n_rows != n_cols)
    {
        // new_b = A^T * b
        MatMultTranspose(A.getRawMatrix(), b.getRawVector(),
                         new_b.getRawVector());
        // new_A = A^T * A
        MatTranspose(A.getRawMatrix(), MAT_INITIAL_MATRIX,
                     &(new_A.getRawMatrix()));
    }
    */
    OGS_FATAL(
        "Normalization operation is not implemented yet for PETSc library! "
        "Program terminated.");
}

References OGS_FATAL.

Referenced by NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::normalizeAandRhs().

matMult() [1/2]

template<typename Matrix , typename Vector >

void MathLib::LinAlg::matMult	(	Matrix const &	A,
		Vector const &	x,
		Vector &	y )

Computes $y = A \cdot x$.

Note

x must not be the same object as y. This restirction has been chosen in order to fulfill the requirements of the respective PETSc function.

Definition at line 123 of file LinAlg.h.

{
    assert(&x != &y);
    y = A * x;
}

matMult() [2/2]

void MathLib::LinAlg::matMult	(	PETScMatrix const &	A,
		PETScVector const &	x,
		PETScVector &	y )

Definition at line 149 of file LinAlg.cpp.

{
    // TODO check sizes
    assert(&x != &y);
    if (!y.getRawVector())
        y.shallowCopy(x);
    MatMult(A.getRawMatrix(), x.getRawVector(), y.getRawVector());
}

References MathLib::PETScMatrix::getRawMatrix(), MathLib::PETScVector::getRawVector(), and MathLib::PETScVector::shallowCopy().

Referenced by NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::calculateNonEquilibriumInitialResiduum(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeResidual(), NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::solve(), and ProcessLib::ComponentTransport::ComponentTransportProcess::solveReactionEquation().

matMultAdd() [1/2]

template<typename Matrix , typename Vector >

void MathLib::LinAlg::matMultAdd	(	Matrix const &	A,
		Vector const &	v1,
		Vector const &	v2,
		Vector &	v3 )

Computes $v_3 = A \cdot v_1 + v_2$.

Note

x must not be the same object as y. This restirction has been chosen in order to fulfill the requirements of the respective PETSc function.

Definition at line 136 of file LinAlg.h.

{
    assert(&v1 != &v3);
    v3 = v2 + A * v1;
}

matMultAdd() [2/2]

void MathLib::LinAlg::matMultAdd	(	PETScMatrix const &	A,
		PETScVector const &	v1,
		PETScVector const &	v2,
		PETScVector &	v3 )

Definition at line 159 of file LinAlg.cpp.

{
    // TODO check sizes
    assert(&v1 != &v3);
    if (!v3.getRawVector())
        v3.shallowCopy(v1);
    MatMultAdd(A.getRawMatrix(), v1.getRawVector(), v2.getRawVector(),
               v3.getRawVector());
}

References MathLib::PETScMatrix::getRawMatrix(), MathLib::PETScVector::getRawVector(), and MathLib::PETScVector::shallowCopy().

Referenced by NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeResidual(), and NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeRhs().

norm()

template<typename MatrixOrVector >

double MathLib::LinAlg::norm	(	MatrixOrVector const &	x,
		MathLib::VecNormType	type )

Definition at line 96 of file LinAlg.h.

{
    switch (type)
    {
        case MathLib::VecNormType::NORM1:
            return norm1(x);
        case MathLib::VecNormType::NORM2:
            return norm2(x);
        case MathLib::VecNormType::INFINITY_N:
            return normMax(x);
        default:
            OGS_FATAL("Invalid norm type given.");
    }
}

References MathLib::INFINITY_N, MathLib::NORM1, norm1(), MathLib::NORM2, norm2(), normMax(), and OGS_FATAL.

Referenced by NumLib::ConvergenceCriterionDeltaX::checkDeltaX(), NumLib::ConvergenceCriterionResidual::checkDeltaX(), NumLib::ConvergenceCriterionResidual::checkResidual(), computeRelativeNorm(), norm1(), norm2(), and normMax().

norm1() [1/2]

template<typename MatrixOrVector >

double MathLib::LinAlg::norm1

(

MatrixOrVector const &

x

)

Computes the Manhattan norm of x.

norm1() [2/2]

template<>

double MathLib::LinAlg::norm1

(

PETScVector const &

x

)

Definition at line 90 of file LinAlg.cpp.

{
    PetscScalar norm = 0.;
    VecNorm(x.getRawVector(), NORM_1, &norm);
    return norm;
}

References MathLib::PETScVector::getRawVector(), and norm().

Referenced by norm().

norm2() [1/2]

template<typename MatrixOrVector >

double MathLib::LinAlg::norm2

(

MatrixOrVector const &

x

)

Computes the Euclidean norm of x.

norm2() [2/2]

template<>

double MathLib::LinAlg::norm2

(

PETScVector const &

x

)

Definition at line 100 of file LinAlg.cpp.

{
    PetscScalar norm = 0.;
    VecNorm(x.getRawVector(), NORM_2, &norm);
    return norm;
}

References MathLib::PETScVector::getRawVector(), and norm().

Referenced by norm().

normMax() [1/2]

template<typename MatrixOrVector >

double MathLib::LinAlg::normMax

(

MatrixOrVector const &

x

)

Computes the maximum norm of x.

normMax() [2/2]

template<>

double MathLib::LinAlg::normMax

(

PETScVector const &

x

)

Definition at line 110 of file LinAlg.cpp.

{
    PetscScalar norm = 0.;
    VecNorm(x.getRawVector(), NORM_INFINITY, &norm);
    return norm;
}

References MathLib::PETScVector::getRawVector(), and norm().

Referenced by norm().

scale() [1/3]

template<typename MatrixOrVector >

void MathLib::LinAlg::scale	(	MatrixOrVector &	x,
		double const	a )

Scales x by a.

Definition at line 51 of file LinAlg.h.

{
    x *= a;
}

scale() [2/3]

void MathLib::LinAlg::scale	(	PETScMatrix &	A,
		PetscScalar const	a )

Definition at line 125 of file LinAlg.cpp.

{
    MatScale(A.getRawMatrix(), a);
}

References MathLib::PETScMatrix::getRawMatrix().

scale() [3/3]

void MathLib::LinAlg::scale	(	PETScVector &	x,
		PetscScalar const	a )

Definition at line 44 of file LinAlg.cpp.

{
    VecScale(x.getRawVector(), a);
}

References MathLib::PETScVector::getRawVector().

Referenced by NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeResidual(), NumLib::TimeDiscretizedODESystem< ODESystemTag::FirstOrderImplicitQuasilinear, NonlinearSolverTag::Newton >::getResidual(), NumLib::BackwardEuler::getWeightedOldX(), and ProcessLib::PhaseField::PhaseFieldProcess< DisplacementDim >::postNonLinearSolverConcreteProcess().

set()

void MathLib::LinAlg::set	(	PETScVector &	x,
		PetscScalar const	a )

Definition at line 32 of file LinAlg.cpp.

{
    VecSet(x.getRawVector(), a);
}

References MathLib::PETScVector::getRawVector().

Referenced by ProcessLib::SmallDeformation::writeMaterialForces().

setLocalAccessibleVector()

void MathLib::LinAlg::setLocalAccessibleVector

(

PETScVector const &

x

)

Set local accessible vector in order to get entries. Call this before call operator[] or get(...) of x.

Definition at line 27 of file LinAlg.cpp.

{
    x.setLocalAccessibleVector();
}

References MathLib::PETScVector::setLocalAccessibleVector().

Referenced by NumLib::LocalLinearLeastSquaresExtrapolator::calculateResidualElement(), ProcessLib::Process::computeSecondaryVariable(), ChemistryLib::PhreeqcIOData::PhreeqcIO::getConcentration(), NumLib::ConvergenceCriterionPerComponentDeltaX::getDampingFactor(), NumLib::ConvergenceCriterionPerComponentResidual::getDampingFactor(), ProcessLib::Process::postIteration(), ProcessLib::TES::TESProcess::postIterationConcreteProcess(), ProcessLib::LIE::HydroMechanics::HydroMechanicsProcess< GlobalDim >::postTimestepConcreteProcess(), ProcessLib::Process::preIteration(), ProcessLib::Process::preOutput(), ProcessLib::Process::preTimestep(), ChemistryLib::PhreeqcIOData::AqueousSolution::print(), ChemistryLib::PhreeqcIOData::anonymous_namespace{PhreeqcIO.cpp}::setAqueousSolution(), ProcessLib::Process::setInitialConditions(), anonymous_namespace{Process.cpp}::setLocalAccessibleVectors(), NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::solve(), ProcessLib::ComponentTransport::ComponentTransportProcess::solveReactionEquation(), NumLib::transformVariableFromGlobalVector(), NumLib::transformVariableFromGlobalVector(), ProcessLib::HMPhaseField::HMPhaseFieldProcess< DisplacementDim >::updateConstraints(), ProcessLib::PhaseField::PhaseFieldProcess< DisplacementDim >::updateConstraints(), and ProcessLib::SmallDeformation::writeMaterialForces().