umfpack_solve.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-2022 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef UMFPACK_SOLVE_HH
#define UMFPACK_SOLVE_HH
 
#include <vector>
#include <memory>
 
#include "linalg/factorization.hh"
 
namespace Kaskade
{
 
// A class for creating and storing symbolic factorizations.
template<class Scalar,class UMFPackIndex=int>
class UMFSymbolic
{
public:
  //
  UMFSymbolic(std::vector<UMFPackIndex> const & Ap,
              std::vector<UMFPackIndex> const & Ai,
              std::vector<Scalar> const & Az,
              int verbose_ = 0);
 
  ~UMFSymbolic(); 
 
  // the opaque symbolic factorization data computed by UMFPack
  void* getMem() { return mem; }
 
  // compressed column storage index arrays
  UMFPackIndex* getAp() { return &lAp[0]; }
  UMFPackIndex* getAi() { return &lAi[0]; }
 
private:
  void* mem;
  std::vector<UMFPackIndex> lAp;
  std::vector<UMFPackIndex> lAi;
};
 
// A class for creating and storing numeric factorizations.
template<class Scalar, class UMFPackIndex=int>
class UMFFactorizationSpace
{
public:
  UMFFactorizationSpace(std::vector<UMFPackIndex> const& Ap,
                        std::vector<UMFPackIndex> const& Ai,
                        std::vector<Scalar> const& Az);
  ~UMFFactorizationSpace();
 
  // the opaque numeric factorization data computed by UMFPack
  void* getMem() { return mem; }
 
  // compressed column storage index arrays
  UMFPackIndex* getAp() { return symbolic.getAp(); }
  UMFPackIndex* getAi() { return symbolic.getAi(); }
 
private:
  UMFSymbolic<Scalar,UMFPackIndex> symbolic;
  void* mem;
};
// ------------------------------------------------------------------------------------------------
 
template <class Scalar, class UMFPackIndex = int>
class UMFFactorization: public Factorization<Scalar>
{
public:
 
  UMFFactorization(UMFPackIndex n_,
                   std::vector<UMFPackIndex> const& ridx,
                   std::vector<UMFPackIndex> const& cidx,
                   std::vector<Scalar> const& values,
                   FactorizationOptions opt)
  : N(ridx.size())
  , n(n_), Az(N)
  , options(opt)
        {
          assert(cidx.size()==N && values.size()==N);
        
          // Transform triplet form into compressed column
    std::vector<UMFPackIndex> Ap(n+1), Ai(N);
          Kaskade::tripletToCompressedColumn(UMFPackIndex(n), UMFPackIndex(n), N,
                                       ridx, cidx, values, Ap, Ai, Az);
          factorization.reset(new UMFFactorizationSpace<Scalar,UMFPackIndex>(Ap,Ai,Az));
        }
 
  UMFFactorization(UMFPackIndex n_,
                   std::unique_ptr<std::vector<UMFPackIndex>> const ridx,
                   std::unique_ptr<std::vector<UMFPackIndex>> const cidx,
                   std::unique_ptr<std::vector<Scalar>> const values,
                   FactorizationOptions opt)
  : N(ridx->size())
  , n(n_), Az(N)
  , options(opt)
  {
    assert(cidx->size()==N && values->size()==N);
    
    // Transform triplet form into compressed column
    std::vector<UMFPackIndex> Ap(n+1), Ai(N);
    Kaskade::tripletToCompressedColumn(UMFPackIndex(n), UMFPackIndex(n), N,
                                       ridx, cidx, values, Ap, Ai, Az);
    
    // Clear memory of data that is not needed anymore.
    ridx.reset();
    cidx.reset();
    values.reset();
    
    factorization.reset(new UMFFactorizationSpace<Scalar,UMFPackIndex>(Ap,Ai,Az));
  }
                                     
 
  void solve(std::vector<Scalar> const& b, std::vector<Scalar>& x, bool transposed = false) const
  {
    assert(b.size() >= n);
    x.resize(n);
    solve_internal(&b[0],&x[0],transposed);
  }
 
  void solve(Scalar const* b, Scalar* x, bool transposed = false) const
  {
    solve_internal(b,x,transposed);
  }
 
  void solve(std::vector<Scalar> const& b, std::vector<Scalar>& x, int nr, bool transposed = false) const
  {
    assert(b.size()>=n*nr);
    x.resize(nr*n);
    for(int i=0; i<nr; ++i)
      solve_internal(&b[i*n],&x[i*n],transposed);
  }
 
  void solve(std::vector<Scalar>& b) const
  {
    assert(b.size() >= n);
    std::vector<Scalar> x(n);
    solve(b,x);
    std::swap(b,x);
  }
 
  void solve(Scalar* b) const
  {
    std::vector<Scalar> x(n);
    solve(b,&x[0]);
    std::copy(x.begin(),x.end(),b);
  }
 
  virtual size_t size() const
  {
    return n;
  }
 
  
private:
 
  virtual void solve_internal(Scalar const* bp, Scalar* xp, bool transposed) const;
 
  size_t const N;   // number of nonzeroes in A
  size_t const n;   // dimension of A
  std::vector<Scalar> Az;
  std::unique_ptr<UMFFactorizationSpace<Scalar,UMFPackIndex>> factorization;
  FactorizationOptions options;
};
 
void umfpack_solve(std::vector<int> const& ridx,
                   std::vector<int> const& cidx,
                   std::vector<double> const& values,
                   std::vector<double> const& b,
                   std::vector<double>& x);
 
}  // namespace Kaskade
#endif