factorization.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef FACTORIZATION_HH
#define FACTORIZATION_HH
 
#include <sstream>
#include <tuple>
#include <variant>
#include <vector>
 
#include <dune/grid/config.h>
 
#include "utilities/enums.hh"
#include "dune/istl/bcrsmatrix.hh"
 
#include "linalg/factorizationOptions.hh"
#include "linalg/factory.hh"
#include "linalg/triplet.hh"
 
namespace Kaskade
{
 
//---------------------------------------------------------------------
 
template <class Scalar>
class Factorization 
{
public:
  typedef Scalar field_type;
 
  using Options = FactorizationOptions;
 
  struct Info
  {
    struct Mumps
    {
      int negativePivots = 0; 
      int zeroPivots = 0; 
    } mumps;
  };
 
  virtual void solve(std::vector<field_type>& b) const
  {
    assert(b.size() >= size());
    solve(&b[0]);
  }
 
  virtual void solve(field_type* b) const = 0;
 
  virtual void solve(std::vector<Scalar> const& b, std::vector<Scalar>& x, bool transposed=false) const
  {
    assert(b.size() >= size());
    if (x.size() < size())
      x.resize(size());
    solve(&b[0], &x[0]);
  }
 
  virtual void solve(Scalar const* b, Scalar* x, bool transposed=false) const = 0;
 
  
  Factorization() = default;
  Factorization(Options options_) : options(options_) {}
  virtual ~Factorization() {}
  void setVerbose(int verbose_) { options.verbosity = verbose_; }
  int getVerbose() const { return options.verbosity; }
 
  Info const& info() const { return info_; }
 
  virtual size_t size() const = 0;
 
protected:
  Options options;
  Info info_;
};
 
// -------------------------------------------------------------------------------------------------
 
template <class Scalar, class Index>
void tripletToCompressedColumn(Index nRows, Index nCols,
                               size_t nNonZeros,
                               std::vector<Index> const& ridx,
                               std::vector<Index> const& cidx,
                               std::vector<Scalar> const& values,
                               std::vector<Index>& Ap, std::vector<Index>& Ai,
                               std::vector<Scalar>& Az);
 
// -------------------------------------------------------------------------------------------------
 
template <class Scalar>
struct Creator<Factorization<Scalar>>
{
  // We can provide a matrix either with int or with long row and column indices.
  // Concrete factorization classes are responsible for conversion, if needed.
  using Matrices = std::variant<MatrixAsTriplet<Scalar,int> const*,
                                MatrixAsTriplet<Scalar,long> const*,
                                MatrixAsTriplet<Scalar,unsigned long> const*>;
  using Argument = std::tuple<Matrices, FactorizationOptions>;
 
  virtual std::unique_ptr<Factorization<Scalar>> create(Argument a) const = 0;
 
  virtual ~Creator() {}
  
  static std::string lookupFailureHint(DirectType direct)
  {
    std::ostringstream hint;
    hint << "Solver key " << direct << " not registered.\n"
         << "Has the corresponding object file been linked in?\n"
         << "Did you specify the right index type (int/long/unsigned long) in your matrix?\n";
    return hint.str();
  }
};
 
 
//---------------------------------------------------------------------
 
template <class Scalar, class Index>
std::unique_ptr<Factorization<Scalar>> getFactorization(DirectType directType,
                                                        MatrixAsTriplet<Scalar,Index> const& A,
                                                        FactorizationOptions options)
{
  assert(A.nrows()==A.ncols());
  if (A.nnz()==0)
    throw SingularMatrixException("No nonzero entries.",__FILE__,__LINE__);
 
  using Argument = typename Creator<Factorization<Scalar>>::Argument;
  using Matrices = typename Creator<Factorization<Scalar>>::Matrices;
 
  Argument args{Matrices(&A),options};
 
  if (directType == DirectType::ANY)
    return Factory<DirectType,Factorization<Scalar>>::createAny(args);
  else
    return Factory<DirectType,Factorization<Scalar>>::create(directType,args);
}
 
template <class Scalar>
std::unique_ptr<Factorization<Scalar>>
getFactorization(DirectType directType,
                 Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>> const& A,
                 FactorizationOptions options)
{
  using Index = typename Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>::size_type;
  return getFactorization(directType,MatrixAsTriplet<Scalar,Index>(A),options);
}
 
template <class Scalar, int n, int m, class Index>
std::unique_ptr<Factorization<Scalar>>
getFactorization(DirectType directType,
                 NumaBCRSMatrix<Dune::FieldMatrix<Scalar,n,m>,Index> const& A,
                 FactorizationOptions options)
{
  return getFactorization(directType,MatrixAsTriplet<Scalar,Index>(A),options);
}
 
}  // namespace Kaskade
 
#endif