mumpsmpi_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-2011 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef MUMPS_SOLVE_HH
#define MUMPS_SOLVE_HH
 
#include <vector>
#include <iostream>
#include <memory>
 
#include "dmumps_c.h"
#define ICNTL(I) icntl[(I)-1]
#define JOB_INIT -1
#define JOB_END -2
#define USE_COMM_WORLD -987654
 
#include "linalg/factorization.hh"
 
namespace Kaskade
{
 
template <class Scalar,class SparseIndexInt=int, class DIL=int>
class MUMPSFactorization: public Factorization<Scalar,SparseIndexInt>
{
public:
 
 
  MUMPSFactorization(SparseIndexInt n_,
                     std::vector<SparseIndexInt> & ridx,
                     std::vector<SparseIndexInt> & cidx,
                     std::vector<Scalar> & values,
                     MatrixProperties property = MatrixProperties::GENERAL,
                     int verb = 0)
  : N(ridx.size()), n(n_)
  {
    assert(cidx.size()==N && values.size()==N);
 
    // Write debug/status output to stdcout if requested.
    id = MPI::COMM_WORLD.Get_rank ( );
    this->setVerbose(verb);
    if ( (this->getVerbose()>=2) && (id == 0) )
    {
      const char *s = "UNKNOWN";
      if (property==MatrixProperties::GENERAL) s = "GENERAL";
      else if (property==MatrixProperties::SYMMETRIC) s = "SYMMETRIC";
      else if (property==MatrixProperties::SYMMETRICSTRUCTURE) s = "SYMMETRICSTRUCTURE";
      else if (property==MatrixProperties::POSITIVEDEFINITE) s = "POSITIVEDEFINITE";
      std::cout << "MUMPS" << " solver, n=" << n << ", nnz=" << N << ", matrix is " << s << std::endl;
    }
 
    init(property);
    analyze(&ridx[0], &cidx[0]);
    factorize(&values[0]);
   };
 
 
  MUMPSFactorization(SparseIndexInt n_,
                     std::unique_ptr<std::vector<SparseIndexInt>> ridx,
                     std::unique_ptr<std::vector<SparseIndexInt>> cidx,
                     std::unique_ptr<std::vector<Scalar>> values,
                     MatrixProperties property = MatrixProperties::GENERAL,
                     int verb = 0)
  : N(ridx->size()), n(n_)
  {
    assert(cidx->size()==N && values->size()==N);
 
    id = MPI::COMM_WORLD.Get_rank ( );
    this->setVerbose(verb);
    if ( (this->getVerbose()>=2) && (id == 0) )
      {
        char *s = "UNKNOWN";
        if (property==MatrixProperties::GENERAL) s = "GENERAL";
        else if (property==MatrixProperties::SYMMETRIC) s = "SYMMETRIC";
        else if (property==MatrixProperties::SYMMETRICSTRUCTURE) s = "SYMMETRICSTRUCTURE";
        else if (property==MatrixProperties::POSITIVEDEFINITE) s = "POSITIVEDEFINITE";
        std::cout << "MUMPS" << " solver, n=" << n << ", nnz=" << N << ", matrix is " << s << std::endl; 
      }
 
    init(property);
    analyze(&ridx[0], &cidx[0]);
    factorize(&values[0]);
  };
                                     
 
  ~MUMPSFactorization()
  {
    mumps_par.job = JOB_END;
    mumps_par.ICNTL(3) = 0;
    dmumps_c(&mumps_par);
  }
 
  void init(MatrixProperties property)
  {
    mumps_par.job = JOB_INIT;
    mumps_par.par = 1;
    switch (property)
    {
      case MatrixProperties::GENERAL:
        mumps_par.sym = 0;
        break;
      case MatrixProperties::SYMMETRICSTRUCTURE:
        mumps_par.sym = 0;
        break;
      case MatrixProperties::SYMMETRIC:
        mumps_par.sym = 2;
        break;
      case MatrixProperties::POSITIVEDEFINITE:
        mumps_par.sym = 1;
        break;
    }
    mumps_par.comm_fortran = USE_COMM_WORLD;
    dmumps_c(&mumps_par);
    if (Factorization<Scalar>::getVerbose()<1)
      mumps_par.ICNTL(3) = 0;
    mumps_par.ICNTL(14) += 5; // 5% additonal workspace
  }
 
  void analyze(SparseIndexInt *irn, SparseIndexInt *jcn)
  {
    mumps_par.job = 1;
    mumps_par.n = n;
    mumps_par.nz = N;
    if ( id == 0 )
      for (int k=0; k<N; k++)
      {
        irn[k]++;
        jcn[k]++;
      };
    MPI_Barrier(MPI_COMM_WORLD);
    mumps_par.irn = irn;
    mumps_par.jcn = jcn;
    dmumps_c(&mumps_par);
  }
 
  void factorize(Scalar *a)
  {
    mumps_par.job = 2;
    mumps_par.a = a;
    dmumps_c(&mumps_par);
  }
 
  void solve(std::vector<Scalar> const& b, std::vector<Scalar>& x, bool transpose=false) const
  {
    assert(b.size()>=n);
    x.resize(n);
    assert(&x != &b);
    x = b;
    mumps_par.job = 3;
    mumps_par.rhs = &x[0];
    dmumps_c(&mumps_par);
  }
 
  void solve(std::vector<Scalar>& b) const
  {
    mumps_par.job = 3;
    mumps_par.rhs = &b[0];
    dmumps_c(&mumps_par);
  }
  
private:
  int id;
  size_t const N; 
  size_t const n;
  mutable DMUMPS_STRUC_C mumps_par;
};
 
}  // namespace Kaskade
#endif