LinAlg.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <cassert>
 
#include "BaseLib/Error.h"
#include "LinAlgEnums.h"
 
#ifdef USE_PETSC
#include <petscsystypes.h>
 
#include "PETSc/PETScVector.h"
#endif
 
namespace MathLib
{

namespace LinAlg
{
 
// Matrix or Vector

template <typename MatrixOrVector>
void copy(MatrixOrVector const& x, MatrixOrVector& y)
{
    y = x;
}

template <typename MatrixOrVector>
void scale(MatrixOrVector& x, double const a)
{
    x *= a;
}

template <typename MatrixOrVector>
void aypx(MatrixOrVector& y, double const a, MatrixOrVector const& x)
{
    y = a * y + x;
}

template <typename MatrixOrVector>
void axpy(MatrixOrVector& y, double const a, MatrixOrVector const& x)
{
    y += a * x;
}

template <typename MatrixOrVector>
void axpby(MatrixOrVector& y, double const a, double const b,
           MatrixOrVector const& x)
{
    y = a * x + b * y;
}

template <typename MatrixOrVector>
void componentwiseDivide(MatrixOrVector& w, MatrixOrVector const& x,
                         MatrixOrVector const& y);

template <typename MatrixOrVector>
double norm1(MatrixOrVector const& x);

template <typename MatrixOrVector>
double norm2(MatrixOrVector const& x);

template <typename MatrixOrVector>
double normMax(MatrixOrVector const& x);
 
template <typename MatrixOrVector>
double norm(MatrixOrVector const& x, MathLib::VecNormType type)
{
    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.");
    }
}
 
template <typename Matrix>
void finalizeAssembly(Matrix& /*A*/);
 
// Matrix and Vector

template <typename Matrix, typename Vector>
void matMult(Matrix const& A, Vector const& x, Vector& y)
{
    assert(&x != &y);
    y = A * x;
}

template <typename Matrix, typename Vector>
void matMultAdd(Matrix const& A, Vector const& v1, Vector const& v2, Vector& v3)
{
    assert(&v1 != &v3);
    v3 = v2 + A * v1;
}
 
}  // namespace LinAlg
}  // namespace MathLib
 
// Global PETScMatrix/PETScVector //////////////////////////////////////////
#ifdef USE_PETSC
 
namespace MathLib
{
 
class PETScMatrix;
class PETScVector;
 
namespace LinAlg
{
 
// Vector

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);
 
// y = a*y + X
void aypx(PETScVector& y, PetscScalar const a, PETScVector const& x);
 
// y = a*x + y
void axpy(PETScVector& y, PetscScalar const a, PETScVector const& x);
 
// y = a*x + b*y
void axpby(PETScVector& y, PetscScalar const a, PetscScalar const b,
           PETScVector const& x);
double dot(PETScVector const& a, PETScVector const& b);
 
// Matrix
 
void copy(PETScMatrix const& A, PETScMatrix& B);
 
// A = a*A
void scale(PETScMatrix& A, PetscScalar const a);
 
// Y = a*Y + X
void aypx(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X);
 
// Y = a*X + Y
void axpy(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X);
 
// Matrix and Vector
 
// v3 = A*v1 + v2
void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y);
 
// y = A*x
void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
                PETScVector const& v2, PETScVector& v3);
 
// new_A = A^T * A
// new_b = A^T * b
void linearSysNormalize(PETScMatrix const& A, PETScMatrix& new_A,
                        PETScVector const& b, PETScVector& new_b);
 
void finalizeAssembly(PETScMatrix& A);
void finalizeAssembly(PETScVector& x);
 
}  // namespace LinAlg
}  // namespace MathLib
 
// Sparse global EigenMatrix/EigenVector
// //////////////////////////////////////////
#else
 
namespace MathLib
{
 
class EigenMatrix;
class EigenVector;
 
namespace LinAlg
{
 
// Vector

void setLocalAccessibleVector(EigenVector const& x);
 
void set(EigenVector& x, double const a);
 
void copy(EigenVector const& x, EigenVector& y);
 
void scale(EigenVector& x, double const a);
 
// y = a*y + X
void aypx(EigenVector& y, double const a, EigenVector const& x);
 
// y = a*x + y
void axpy(EigenVector& y, double const a, EigenVector const& x);
 
// y = a*x + b*y
void axpby(EigenVector& y, double const a, double const b,
           EigenVector const& x);
double dot(EigenVector const& a, EigenVector const& b);
 
// Matrix
 
void copy(EigenMatrix const& A, EigenMatrix& B);
 
// A = a*A
void scale(EigenMatrix& A, double const a);
 
// Y = a*Y + X
void aypx(EigenMatrix& Y, double const a, EigenMatrix const& X);
 
// Y = a*X + Y
void axpy(EigenMatrix& Y, double const a, EigenMatrix const& X);
 
// Matrix and Vector
 
// y = A*x
void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y);
 
// v3 = A*v1 + v2
void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
                EigenVector const& v2, EigenVector& v3);
// new_A = A^T * A
// new_b = A^T * b
void linearSysNormalize(EigenMatrix const& A, EigenMatrix& new_A,
                        EigenVector const& b, EigenVector& new_b);
void finalizeAssembly(EigenMatrix& x);
void finalizeAssembly(EigenVector& A);
 
}  // namespace LinAlg
 
}  // namespace MathLib
 
#endif
 
namespace MathLib::LinAlg
{

template <typename VectorType>
double computeRelativeNorm(VectorType const& x, VectorType const& y,
                           MathLib::VecNormType norm_type)
{
    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();
}
}  // namespace MathLib::LinAlg