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)
 
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 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]

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

72{
73 y = a*x + b*y;
74}

◆ 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.

66{
67 // TODO check sizes
68 VecAXPBY(y.getRawVector(), a, b, x.getRawVector());
69}
PETSc_Vec & getRawVector()
Exposes the underlying PETSc vector.

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

◆ axpy() [1/3]

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

Computes \( y = a \cdot x + y \).

Definition at line 63 of file LinAlg.h.

64{
65 y += a*x;
66}

◆ axpy() [2/3]

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

Definition at line 139 of file LinAlg.cpp.

140{
141 // TODO check sizes
142 // TODO sparsity pattern, currently they are assumed to be different (slow)
143 MatAXPY(Y.getRawMatrix(), a, X.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
144}
Mat & getRawMatrix()
Get matrix reference.
Definition PETScMatrix.h:86

References MathLib::PETScMatrix::getRawMatrix(), and set().

◆ axpy() [3/3]

◆ aypx() [1/3]

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

Computes \( y = a \cdot y + x \).

Definition at line 56 of file LinAlg.h.

57{
58 y = a*y + x;
59}

◆ aypx() [2/3]

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

Definition at line 131 of file LinAlg.cpp.

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

References MathLib::PETScMatrix::getRawMatrix(), and set().

◆ aypx() [3/3]

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

Definition at line 50 of file LinAlg.cpp.

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

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

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

◆ componentwiseDivide() [1/2]

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.

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

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

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

288{
289 if (norm_type == MathLib::VecNormType::INVALID)
290 {
291 OGS_FATAL("An invalid norm type has been passed");
292 }
293
294 // Stores \f$diff = x-y\f$.
295 VectorType diff;
296 MathLib::LinAlg::copy(x, diff);
297 MathLib::LinAlg::axpy(diff, -1.0, y);
298
299 const double norm_diff = MathLib::LinAlg::norm(diff, norm_type);
300
301 const double norm_x = MathLib::LinAlg::norm(x, norm_type);
302 if (norm_x > std::numeric_limits<double>::epsilon())
303 {
304 return norm_diff / norm_x;
305 }
306
307 // Both of norm_x and norm_diff are close to zero
308 if (norm_diff < std::numeric_limits<double>::epsilon())
309 {
310 return 1.0;
311 }
312
313 // Only norm_x is close to zero
314 return norm_diff / std::numeric_limits<double>::epsilon();
315}
#define OGS_FATAL(...)
Definition Error.h:26
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition LinAlg.h:94
void copy(PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:37
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition LinAlg.cpp:57

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

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

◆ copy() [1/3]

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

Copies x to y.

Definition at line 42 of file LinAlg.h.

43{
44 y = x;
45}

◆ copy() [2/3]

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

Definition at line 119 of file LinAlg.cpp.

120{
121 B = A;
122}

◆ copy() [3/3]

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

176{
178}
void finalizeAssembly()
Perform MPI collection of assembled entries in buffer.

References MathLib::PETScVector::finalizeAssembly().

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

121{
122 assert(&x != &y);
123 y = A*x;
124}

References set().

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

134{
135 assert(&v1 != &v3);
136 v3 = v2 + A*v1;
137}

References set().

◆ 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.

161{
162 // TODO check sizes
163 assert(&v1 != &v3);
164 if (!v3.getRawVector())
165 v3.shallowCopy(v1);
166 MatMultAdd(A.getRawMatrix(), v1.getRawVector(), v2.getRawVector(),
167 v3.getRawVector());
168}

References MathLib::PETScMatrix::getRawMatrix(), and set().

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

◆ norm()

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

Definition at line 94 of file LinAlg.h.

95{
96 switch (type) {
98 return norm1(x);
100 return norm2(x);
102 return normMax(x);
103 default:
104 OGS_FATAL("Invalid norm type given.");
105 }
106}
double norm1(PETScVector const &x)
Definition LinAlg.cpp:90
double normMax(PETScVector const &x)
Definition LinAlg.cpp:110
double norm2(PETScVector const &x)
Definition LinAlg.cpp:100

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

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

◆ norm1() [1/2]

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.

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

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

Referenced by norm().

◆ norm2() [1/2]

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.

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

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

Referenced by norm().

◆ normMax() [1/2]

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.

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

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

Referenced by norm().

◆ scale() [1/3]

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

Scales x by a.

Definition at line 49 of file LinAlg.h.

50{
51 x *= a;
52}

◆ scale() [2/3]

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

Definition at line 125 of file LinAlg.cpp.

126{
127 MatScale(A.getRawMatrix(), a);
128}

References MathLib::PETScMatrix::getRawMatrix(), and set().

◆ scale() [3/3]

◆ set()

◆ setLocalAccessibleVector()

void MathLib::LinAlg::setLocalAccessibleVector ( PETScVector const & x)