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 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 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().

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

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

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

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.

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

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

299{
300 if (norm_type == MathLib::VecNormType::INVALID)
301 {
302 OGS_FATAL("An invalid norm type has been passed");
303 }
304
305 // Stores \f$diff = x-y\f$.
306 VectorType diff;
307 MathLib::LinAlg::copy(x, diff);
308 MathLib::LinAlg::axpy(diff, -1.0, y);
309
310 const double norm_diff = MathLib::LinAlg::norm(diff, norm_type);
311
312 const double norm_x = MathLib::LinAlg::norm(x, norm_type);
313 if (norm_x > std::numeric_limits<double>::epsilon())
314 {
315 return norm_diff / norm_x;
316 }
317
318 // Both of norm_x and norm_diff are close to zero
319 if (norm_diff < std::numeric_limits<double>::epsilon())
320 {
321 return 1.0;
322 }
323
324 // Only norm_x is close to zero
325 return norm_diff / std::numeric_limits<double>::epsilon();
326}
#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(), 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 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]

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

Definition at line 37 of file LinAlg.cpp.

38{
39 if (!y.getRawVector())
40 y.shallowCopy(x);
41 VecCopy(x.getRawVector(), y.getRawVector());
42}
void shallowCopy(const PETScVector &v)

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::PhaseField::PhaseFieldProcess< DisplacementDim >::updateConstraints(), and NumLib::StaggeredCoupling::updatePreviousSolution().

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

204{
206}
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 170 of file LinAlg.cpp.

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

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

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

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

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

◆ 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(), 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 94 of file LinAlg.h.

95{
96 switch (type)
97 {
99 return norm1(x);
101 return norm2(x);
103 return normMax(x);
104 default:
105 OGS_FATAL("Invalid norm type given.");
106 }
107}
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, 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.

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

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.

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

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.

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

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 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().

◆ scale() [3/3]

◆ set()

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

Definition at line 32 of file LinAlg.cpp.

33{
34 VecSet(x.getRawVector(), a);
35}

References MathLib::PETScVector::getRawVector().

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

◆ setLocalAccessibleVector()