umfpack_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_UMFPACK_SOLVER_H_
#define __HERMES_COMMON_UMFPACK_SOLVER_H_
#include "config.h"
#ifdef WITH_UMFPACK
#include "linear_matrix_solver.h"
#include "matrix.h"

using namespace Hermes::Algebra;

namespace Hermes
{
  namespace Solvers
  {
    template <typename Scalar> class HERMES_API UMFPackLinearMatrixSolver;
    template <typename Scalar> class HERMES_API UMFPackIterator;
  }

  namespace Algebra
  {
    using namespace Hermes::Solvers;
    template <typename Scalar>
    class HERMES_API CSCMatrix : public SparseMatrix<Scalar>
    {
    public:
      void create(unsigned int size, unsigned int nnz, int* ap, int* ai, Scalar* ax);

      CSCMatrix();
      CSCMatrix(unsigned int size);
      virtual ~CSCMatrix();
      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_matrix(CSCMatrix<Scalar>* mat);
      virtual void add_to_diagonal_blocks(int num_stages, CSCMatrix<Scalar>* 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, CSCMatrix<Scalar>* mat);
      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;

      // Applies the matrix to vector_in and saves result to vector_out.
      void multiply_with_vector(Scalar* vector_in, Scalar* vector_out);
      // Multiplies matrix with a Scalar.
      void multiply_with_Scalar(Scalar value);

      // Duplicates a matrix (including allocation).
      CSCMatrix* duplicate();
      // Exposes pointers to the CSC arrays.
      int *get_Ap();
      // Exposes pointers to the CSC arrays.
      int *get_Ai();
      // Exposes pointers to the CSC arrays.
      Scalar *get_Ax();

    protected:
      // UMFPack specific data structures for storing the system matrix (CSC format).
      Scalar *Ax;
      int *Ai;
      int *Ap;
      unsigned int nnz;
      template <typename T> friend class Hermes::Solvers::UMFPackLinearMatrixSolver;
      template <typename T> friend class Hermes::Solvers::UMFPackIterator;
      template<typename T> friend SparseMatrix<T>*  create_matrix();
    };

    template <typename Scalar>
    class HERMES_API UMFPackMatrix : public CSCMatrix<Scalar>
    {
      template <typename T> friend class Hermes::Solvers::UMFPackLinearMatrixSolver;
      template <typename T> friend class Hermes::Solvers::UMFPackIterator;
      template<typename T> friend SparseMatrix<T>*  create_matrix();
    };

    template <typename Scalar>
    class HERMES_API UMFPackVector : public Vector<Scalar>
    {
    public:
      UMFPackVector();
      UMFPackVector(unsigned int size);
      virtual ~UMFPackVector();
      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");

      Scalar *get_c_array();

    protected:
      Scalar *v;
      template <typename T> friend class Hermes::Solvers::UMFPackLinearMatrixSolver;
      template <typename T> friend class Hermes::Solvers::UMFPackIterator;
      template<typename T> friend Vector<T>* Hermes::Algebra::create_vector();
    };
  }
  namespace Solvers
  {
    template <typename Scalar>
    class HERMES_API UMFPackLinearMatrixSolver : public DirectSolver<Scalar>
    {
    public:
      UMFPackLinearMatrixSolver(UMFPackMatrix<Scalar> *m, UMFPackVector<Scalar> *rhs);
      virtual ~UMFPackLinearMatrixSolver();
      virtual bool solve();
      virtual int get_matrix_size();

      UMFPackMatrix<Scalar> *m;
      UMFPackVector<Scalar> *rhs;

      void *symbolic;
      void *numeric;  

      void free_factorization_data();
      bool setup_factorization();
      template <typename T> friend class Hermes::Algebra::CSCMatrix;
      template <typename T> friend class Hermes::Algebra::UMFPackMatrix;
      template <typename T> friend class Hermes::Algebra::UMFPackVector;
      template<typename T> friend LinearMatrixSolver<T>* create_linear_solver(Matrix<T>* matrix, Vector<T>* rhs);
      void check_status(const char *fn_name, int status);
    };

    template <typename Scalar>
    class UMFPackIterator
    {
    protected:
      UMFPackIterator(CSCMatrix<Scalar>* mat);
      bool init();
      void get_current_position(int& i, int& j, Scalar& val);
      bool move_to_position(int i, int j);
      bool move_ptr();
      void add_to_current_position(Scalar val);
    protected:
      int size;
      int nnz;
      int* Ai;
      int* Ap;
      Scalar* Ax;
      int Ai_pos;
      int Ap_pos;
      template <typename T> friend class Hermes::Algebra::CSCMatrix;
      template <typename T> friend class Hermes::Algebra::UMFPackMatrix;
      template <typename T> friend class Hermes::Algebra::UMFPackVector;
    };
  }
}
#endif
#endif