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://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 "linear_matrix_solver.h"
#include "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 SparseMatrix<Scalar>
    {
    public:
      MumpsMatrix();
      virtual ~MumpsMatrix();

      virtual void alloc();
      virtual void free();
      virtual Scalar get(unsigned int m, unsigned int n);
      virtual void zero();
      virtual void add(unsigned int m, unsigned int n, Scalar v);
      virtual void add_to_diagonal(Scalar v);
      virtual void add(unsigned int m, unsigned int n, Scalar **mat, int *rows, int *cols);
      virtual bool dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt = DF_MATLAB_SPARSE, char* number_format = "%lf");
      virtual unsigned int get_matrix_size() const;
      virtual unsigned int get_nnz() const;
      virtual double get_fill_in() const;
      virtual void add_matrix(MumpsMatrix* mat);
      virtual void add_to_diagonal_blocks(int num_stages, MumpsMatrix* mat);
      virtual void add_sparse_to_diagonal_blocks(int num_stages, SparseMatrix<Scalar>* mat);
      virtual void add_as_block(unsigned int i, unsigned int j, MumpsMatrix* mat);

      void multiply_with_vector(Scalar* vector_in, Scalar* vector_out);
      void multiply_with_Scalar(Scalar value);
      void create(unsigned int size, unsigned int nnz, int* ap, int* ai, Scalar* ax);
      MumpsMatrix* duplicate();

    protected:
      unsigned int nnz;          
      int *irn;         
      int *jcn;         
      typename mumps_type<Scalar>::mumps_Scalar *Ax; 
      int *Ai;          
      unsigned int *Ap;          

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

    template <typename Scalar>
    class MumpsVector : public Vector<Scalar>
    {
    public:
      MumpsVector();
      virtual ~MumpsVector();

      virtual void alloc(unsigned int ndofs);
      virtual void free();
      virtual Scalar get(unsigned int idx);
      virtual void extract(Scalar *v) const;
      virtual void zero();
      virtual void change_sign();
      virtual void set(unsigned int idx, Scalar y);
      virtual void add(unsigned int idx, Scalar y);
      virtual void add(unsigned int n, unsigned int *idx, Scalar *y);
      virtual void add_vector(Vector<Scalar>* vec);
      virtual void add_vector(Scalar* vec);
      virtual bool dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt = DF_MATLAB_SPARSE, char* number_format = "%lf");

    protected:
      // MUMPS specific data structures for storing the rhs.
      Scalar *v;

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

      virtual bool solve();
      virtual int get_matrix_size();

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

      bool setup_factorization();

      typename mumps_type<Scalar>::mumps_struct param;

      bool check_status();

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