OGS
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)
double dot (PETScVector const &a, PETScVector const &b)
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 65 of file LinAlg.h.

67{
68 y = a * x + b * y;
69}

◆ axpby() [2/2]

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

Definition at line 57 of file LinAlg.cpp.

59{
60 // TODO check sizes
61 VecAXPBY(y.getRawVector(), a, b, x.getRawVector());
62}
PETSc_Vec & getRawVector()
Exposes the underlying PETSc vector.

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 58 of file LinAlg.h.

59{
60 y += a * x;
61}

◆ axpy() [2/3]

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

Definition at line 141 of file LinAlg.cpp.

142{
143 auto const pattern = nonzeroPatternStructure(X.getSparsityColumnIndices(),
145 MatAXPY(Y.getRawMatrix(), a, X.getRawMatrix(), pattern);
146}
Mat & getRawMatrix()
Get matrix reference.
Definition PETScMatrix.h:85
std::vector< PetscInt > const & getSparsityColumnIndices() const
Column indices cached at preallocation (empty on the fallback path).
Definition PETScMatrix.h:79
MatStructure nonzeroPatternStructure(std::vector< PetscInt > const &col_a, std::vector< PetscInt > const &col_b)
double const GaussLegendre< 1 >::X[1]

References MathLib::GaussLegendre< 1 >::X, MathLib::PETScMatrix::getRawMatrix(), MathLib::PETScMatrix::getSparsityColumnIndices(), and MathLib::nonzeroPatternStructure().

◆ axpy() [3/3]

◆ 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 51 of file LinAlg.h.

52{
53 y = a * y + x;
54}

◆ aypx() [2/3]

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

Definition at line 133 of file LinAlg.cpp.

134{
135 auto const pattern = nonzeroPatternStructure(X.getSparsityColumnIndices(),
137 MatAYPX(Y.getRawMatrix(), a, X.getRawMatrix(), pattern);
138}

References MathLib::GaussLegendre< 1 >::X, MathLib::PETScMatrix::getRawMatrix(), MathLib::PETScMatrix::getSparsityColumnIndices(), and MathLib::nonzeroPatternStructure().

◆ aypx() [3/3]

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

Definition at line 43 of file LinAlg.cpp.

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

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 83 of file LinAlg.cpp.

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

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
xVector x
yVector y
norm_typeThe norm type of global vector
Returns
\( \|x-y\|/\|x\| \).

Definition at line 298 of file LinAlg.h.

300{
301 if (norm_type == MathLib::VecNormType::INVALID)
302 {
303 OGS_FATAL("An invalid norm type has been passed");
304 }
305
306 // Stores \f$diff = x-y\f$.
307 VectorType diff;
308 MathLib::LinAlg::copy(x, diff);
309 MathLib::LinAlg::axpy(diff, -1.0, y);
310
311 const double norm_diff = MathLib::LinAlg::norm(diff, norm_type);
312
313 const double norm_x = MathLib::LinAlg::norm(x, norm_type);
314 if (norm_x > std::numeric_limits<double>::epsilon())
315 {
316 return norm_diff / norm_x;
317 }
318
319 // Both of norm_x and norm_diff are close to zero
320 if (norm_diff < std::numeric_limits<double>::epsilon())
321 {
322 return 1.0;
323 }
324
325 // Only norm_x is close to zero
326 return norm_diff / std::numeric_limits<double>::epsilon();
327}
#define OGS_FATAL(...)
Definition Error.h:10
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition LinAlg.h:89
void copy(PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:30
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition LinAlg.cpp:50

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 37 of file LinAlg.h.

38{
39 y = x;
40}

◆ copy() [2/3]

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

Definition at line 121 of file LinAlg.cpp.

122{
123 B = A;
124}

◆ copy() [3/3]

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

Definition at line 30 of file LinAlg.cpp.

31{
32 if (!y.getRawVector())
33 y.shallowCopy(x);
34 VecCopy(x.getRawVector(), y.getRawVector());
35}
void shallowCopy(const PETScVector &v)

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

Referenced by NumLib::AndersonAcceleration::accelerate(), ProcessLib::AssemblyMixin< Process >::assembleOnBulkMesh(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeA(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeJacobian(), computeRelativeNorm(), NumLib::MatrixTranslatorGeneral< ODESystemTag::FirstOrderImplicitQuasilinear >::computeResidual(), ProcessLib::computeResiduum(), 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().

◆ dot()

double MathLib::LinAlg::dot ( PETScVector const & a,
PETScVector const & b )

Computes the dot product \( a \cdot b \) over the global vectors, i.e. ghost entries are not counted twice.

Note
Collective operation; all ranks obtain the same result.

Definition at line 64 of file LinAlg.cpp.

65{
66 PetscScalar result;
67 VecDot(a.getRawVector(), b.getRawVector(), &result);
68 // PetscRealPart is a no-op for real builds and keeps this compiling for a
69 // complex-scalar PETSc build (where PetscScalar is std::complex).
70 return static_cast<double>(PetscRealPart(result));
71}

References MathLib::PETScVector::getRawVector().

Referenced by NumLib::AndersonAcceleration::accelerate().

◆ finalizeAssembly() [1/3]

template<typename Matrix>
void MathLib::LinAlg::finalizeAssembly ( Matrix & )

◆ finalizeAssembly() [2/3]

◆ finalizeAssembly() [3/3]

void MathLib::LinAlg::finalizeAssembly ( PETScVector & x)

Definition at line 205 of file LinAlg.cpp.

206{
208}
void finalizeAssembly()
Perform MPI collection of assembled entries in buffer.

References MathLib::PETScVector::finalizeAssembly().

◆ linearSysNormalize()

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

Definition at line 172 of file LinAlg.cpp.

174{
175 // The following block is deactivated, because there is no tests yet for the
176 // normalization operation in PETSc. This will be a task for later.
177 /*
178 assert(&A != &new_A);
179 assert(&b != &new_b);
180
181 PetscInt n_rows(0);
182 PetscInt n_cols(0);
183 MatGetSize(A.getRawMatrix(), &n_rows, &n_cols);
184 // only when A matrix is not square
185 if (n_rows != n_cols)
186 {
187 // new_b = A^T * b
188 MatMultTranspose(A.getRawMatrix(), b.getRawVector(),
189 new_b.getRawVector());
190 // new_A = A^T * A
191 MatTranspose(A.getRawMatrix(), MAT_INITIAL_MATRIX,
192 &(new_A.getRawMatrix()));
193 }
194 */
195 OGS_FATAL(
196 "Normalization operation is not implemented yet for PETSc library! "
197 "Program terminated.");
198}

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 116 of file LinAlg.h.

117{
118 assert(&x != &y);
119 y = A * x;
120}

◆ matMult() [2/2]

◆ 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 129 of file LinAlg.h.

130{
131 assert(&v1 != &v3);
132 v3 = v2 + A * v1;
133}

◆ matMultAdd() [2/2]

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

◆ norm()

template<typename MatrixOrVector>
double MathLib::LinAlg::norm ( MatrixOrVector const & x,
MathLib::VecNormType type )

Definition at line 89 of file LinAlg.h.

90{
91 switch (type)
92 {
94 return norm1(x);
96 return norm2(x);
98 return normMax(x);
99 default:
100 OGS_FATAL("Invalid norm type given.");
101 }
102}
double norm1(PETScVector const &x)
Definition LinAlg.cpp:92
double normMax(PETScVector const &x)
Definition LinAlg.cpp:112
double norm2(PETScVector const &x)
Definition LinAlg.cpp:102

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 92 of file LinAlg.cpp.

93{
94 PetscScalar norm = 0.;
95 VecNorm(x.getRawVector(), NORM_1, &norm);
96 return norm;
97}

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 102 of file LinAlg.cpp.

103{
104 PetscScalar norm = 0.;
105 VecNorm(x.getRawVector(), NORM_2, &norm);
106 return norm;
107}

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 112 of file LinAlg.cpp.

113{
114 PetscScalar norm = 0.;
115 VecNorm(x.getRawVector(), NORM_INFINITY, &norm);
116 return norm;
117}

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 44 of file LinAlg.h.

45{
46 x *= a;
47}

◆ scale() [2/3]

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

Definition at line 127 of file LinAlg.cpp.

128{
129 MatScale(A.getRawMatrix(), a);
130}

References MathLib::PETScMatrix::getRawMatrix().

◆ scale() [3/3]

◆ set()

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

Definition at line 25 of file LinAlg.cpp.

26{
27 VecSet(x.getRawVector(), a);
28}

References MathLib::PETScVector::getRawVector().

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

◆ setLocalAccessibleVector()