mumps_solver.h

Go to the documentation of this file.

// This file is part of HermesCommon
//
// Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
// Email: hpfem-group@unr.edu, home page: http://www.hpfem.org/.
//
// Hermes2D is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation; either version 2 of the License,
// or (at your option) any later version.
//
// Hermes2D is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Hermes2D; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#ifndef __HERMES_COMMON_MUMPS_SOLVER_H_
#define __HERMES_COMMON_MUMPS_SOLVER_H_
#include "config.h"
#ifdef WITH_MUMPS
#include "solvers/linear_matrix_solver.h"
#include "algebra/matrix.h"

extern "C"
{
#include <mumps_c_types.h>
#include <dmumps_c.h>
#include <zmumps_c.h>
}

#ifdef WITH_MPI
#include <mpi.h>
#endif

namespace Hermes
{
  namespace Solvers
  {
    template <typename Scalar> class MumpsSolver;
  }
}
namespace Hermes
{
  namespace Algebra
  {
    template <typename Scalar> struct mumps_type;

    template <>
    struct mumps_type < std::complex<double> >
    {
      typedef ZMUMPS_STRUC_C mumps_struct;
      typedef ZMUMPS_COMPLEX mumps_Scalar;
    };

    template <>
    struct mumps_type < double >
    {
      typedef DMUMPS_STRUC_C mumps_struct;
      typedef double mumps_Scalar;
    };

    template <typename Scalar>
    class MumpsMatrix : public CSCMatrix < Scalar >
    {
    public:
      MumpsMatrix();
      virtual ~MumpsMatrix();

      void alloc_data();
      void free();
      Scalar get(unsigned int m, unsigned int n) const;
      void zero();

      void add(unsigned int m, unsigned int n, Scalar v);

      void export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char* number_format = "%lf");

      void import_from_file(const char *filename, const char *var_name, MatrixExportFormat fmt);

      void multiply_with_Scalar(Scalar value);

      void multiply_with_vector(Scalar* vector_in, Scalar*& vector_out, bool vector_out_initialized = false) const;

      void create(unsigned int size, unsigned int nnz, int* ap, int* ai, Scalar* ax);

      CSMatrix<Scalar>* duplicate() const;

    protected:
      int *irn;
      int *jcn;
      typename mumps_type<Scalar>::mumps_Scalar *Ax;

      friend class Solvers::MumpsSolver < Scalar > ;
      template<typename T> friend SparseMatrix<T>*  create_matrix();
    };
  }

  namespace Solvers
  {
    template <typename Scalar>
    class HERMES_API MumpsSolver : public DirectSolver < Scalar >
    {
    public:
      MumpsSolver(MumpsMatrix<Scalar> *m, SimpleVector<Scalar> *rhs);
      virtual ~MumpsSolver();
      void free();

      virtual void solve();
      virtual int get_matrix_size();

      MumpsMatrix<Scalar> *m;
      SimpleVector<Scalar> *rhs;

      bool setup_factorization();

      typename mumps_type<Scalar>::mumps_struct param;

      bool check_status();

      bool reinit();
    private:
      //wrapper around dmums_c or zmumps_c
      void mumps_c(typename mumps_type<Scalar>::mumps_struct * param);

      bool inited;

      int icntl_14;
      static const int init_icntl_14 = 1;
      static const int max_icntl_14 = 128;
    };
  }
}
#endif
#endif