HypreParMatrix.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

/// This HYPRE library interface has been taken originally from MFEM and modified
/// to suit the needs for the MOAB library.
/// Modified by: Vijay Mahadevan

#ifndef MOAB_HYPREPARMATRIX
#define MOAB_HYPREPARMATRIX

#include "moab/MOABConfig.h"
#include "moab/Core.hpp"

#ifdef MOAB_HAVE_EIGEN3
#include 
#include 
#else
#error Configure with Eigen3 enabled
#endif

#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"

// hypre header files
// #include "HYPRE.h"
// #include "HYPRE_parcsr_ls.h"
// #include "HYPRE_MatvecFunctions.h"

#include "seq_mv.h"
#include "HypreParVector.hpp"
#include "_hypre_parcsr_mv.h"
#include "_hypre_parcsr_ls.h"
#include "interpreter.h"
#include "HYPRE_MatvecFunctions.h"

#include "temp_multivector.h"

#ifdef HYPRE_COMPLEX
#error "MOAB does not work with HYPRE's complex numbers support"
#endif

#include "hypre_parcsr.hpp"

namespace moab
{

namespace internal
{

    // Convert a HYPRE_Int to int
    inline int to_int( HYPRE_Int i )
    {
#ifdef HYPRE_BIGINT
        MFEM_ASSERT( HYPRE_Int( int( i ) ) == i, "overflow converting HYPRE_Int to int" );
#endif
        return int( i );
    }

}  // namespace internal

/// Wrapper for hypre's ParCSR matrix class
class HypreParMatrix
{
  private:
    /// The actual object
    HYPRE_IJMatrix A;
    hypre_ParCSRMatrix* A_parcsr;

    mv_InterfaceInterpreter* interpreter;

    // Does the object own the pointer A?
    char ParCSROwner;

    // Store the dimensions of the matrix
    int height, gnrows;
    int width, gncols;

    moab::ParallelComm* pcomm;
    char initialized;

    // Initialize with defaults. Does not initialize inherited members.
    void Init();

    // Delete all owned data. Does not perform re-initialization with defaults.
    void Destroy();

    friend class HypreSolver;

  public:
    typedef Eigen::VectorXd Vector;
    typedef Eigen::SparseMatrix< double, Eigen::RowMajor > MOABSparseMatrix;
    // typedef template Eigen::ArrayXX Array2D;

  public:
    /// An empty matrix to be used as a reference to an existing matrix
    HypreParMatrix( moab::ParallelComm* p_comm );

    /// Converts hypre's format to HypreParMatrix
    HypreParMatrix( HYPRE_IJMatrix a )
    {
        Init();
        A      = a;
        height = GetNumRows();
        width  = GetNumCols();
    }

    /** Creates block-diagonal square parallel matrix. Diagonal is given by diag
        which must be in CSR format (finalized). The new HypreParMatrix does not
        take ownership of any of the input arrays. */
    HypreParMatrix( moab::ParallelComm* p_comm, HYPRE_Int glob_size, HYPRE_Int* row_starts, HYPRE_Int nnz_pr_diag = 0 );

    /** Creates block-diagonal rectangular parallel matrix. Diagonal is given by
        diag which must be in CSR format (finalized). The new HypreParMatrix does
        not take ownership of any of the input arrays. */
    HypreParMatrix( moab::ParallelComm* p_comm,
                    HYPRE_Int global_num_rows,
                    HYPRE_Int global_num_cols,
                    HYPRE_Int* row_starts,
                    HYPRE_Int* col_starts,
                    HYPRE_Int nnz_pr_diag    = 0,
                    HYPRE_Int onz_pr_diag    = 0,
                    HYPRE_Int nnz_pr_offdiag = 0 );

    /** Creates a general parallel matrix from a local CSR matrix on each
        processor described by the I, J and data arrays. The local matrix should
        be of size (local) nrows by (global) glob_ncols. The new parallel matrix
        contains copies of all input arrays (so they can be deleted). */
    // HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_Int glob_nrows,
    //                HYPRE_Int glob_ncols, int *I, HYPRE_Int *J,
    //                double *data, HYPRE_Int *rows, HYPRE_Int *cols);

    /// Make this HypreParMatrix a reference to 'master'
    void MakeRef( const HypreParMatrix& master );

    /// MPI communicator
    moab::ParallelComm* GetParallelCommunicator() const
    {
        return pcomm;
    }

    /// Typecasting to hypre's HYPRE_IJMatrix*
    operator HYPRE_IJMatrix()
    {
        return A;
    }
    /// Typecasting to hypre's HYPRE_ParCSRMatrix, a.k.a. void *
    operator HYPRE_ParCSRMatrix()
    {
        return (HYPRE_ParCSRMatrix)A_parcsr;
    }

    /// Changes the ownership of the matrix
    HYPRE_IJMatrix StealData();
    /// Changes the ownership of the matrix to A
    void StealData( HypreParMatrix& A );

    /** If the HypreParMatrix does not own the row-starts array, make a copy of
        it that the HypreParMatrix will own. If the col-starts array is the same
        as the row-starts array, col-starts is also replaced. */
    // void CopyRowStarts();
    /** If the HypreParMatrix does not own the col-starts array, make a copy of
        it that the HypreParMatrix will own. If the row-starts array is the same
        as the col-starts array, row-starts is also replaced. */
    // void CopyColStarts();

    /// Returns the global number of nonzeros
    inline HYPRE_Int NNZ()
    {
        return A_parcsr->num_nonzeros;
    }
    /// Returns the row partitioning
    inline HYPRE_Int* RowPart()
    {
        return A->row_partitioning;
    }
    /// Returns the column partitioning
    inline HYPRE_Int* ColPart()
    {
        return A->col_partitioning;
    }
    /// Returns the global number of rows
    inline HYPRE_Int M()
    {
        return A->global_num_rows;
    }
    /// Returns the global number of columns
    inline HYPRE_Int N()
    {
        return A->global_num_cols;
    }
    /// Returns the global number of rows
    inline HYPRE_Int FirstRow()
    {
        return A->global_first_row;
    }
    /// Returns the global number of columns
    inline HYPRE_Int FirstCol()
    {
        return A->global_first_col;
    }

    /// Get the local diagonal of the matrix.
    void GetDiag( Vector& diag ) const;
    /// Get the local diagonal block. NOTE: 'diag' will not own any data.
    void GetDiag( MOABSparseMatrix& diag ) const;
    /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
    void GetOffd( MOABSparseMatrix& offd, HYPRE_Int*& cmap ) const;

    /** Split the matrix into M x N equally sized blocks of parallel matrices.
        The size of 'blocks' must already be set to M x N. */
    // void GetBlocks(Eigen::Array &blocks,
    //                bool interleaved_rows = false,
    //                bool interleaved_cols = false) const;

    /// Returns the transpose of *this
    // HypreParMatrix * Transpose();

    /** Destroy and resize block-diagonal square parallel matrix. Diagonal is given by diag
        which must be in CSR format (finalized). The new HypreParMatrix does not
        take ownership of any of the input arrays. */
    void resize( HYPRE_Int glob_size, HYPRE_Int* row_starts, HYPRE_Int nnz_pr_diag = 0, HYPRE_Int nnz_pr_offdiag = 0 );

    /** Destroy and resize block-diagonal rectangular parallel matrix. Diagonal is given by
        diag which must be in CSR format (finalized). The new HypreParMatrix does
        not take ownership of any of the input arrays. */
    void resize( HYPRE_Int global_num_rows,
                 HYPRE_Int global_num_cols,
                 HYPRE_Int* row_starts,
                 HYPRE_Int* col_starts,
                 HYPRE_Int* nnz_pr_diag   = NULL,
                 HYPRE_Int* onz_pr_diag   = NULL,
                 HYPRE_Int nnz_pr_offdiag = 0 );

    /// Returns the number of rows in the diagonal block of the ParCSRMatrix
    int GetNumRows() const
    {
        return internal::to_int( hypre_CSRMatrixNumRows( hypre_ParCSRMatrixDiag( A_parcsr ) ) );
    }

    /// Returns the number of columns in the diagonal block of the ParCSRMatrix
    int GetNumCols() const
    {
        return internal::to_int( hypre_CSRMatrixNumCols( hypre_ParCSRMatrixDiag( A_parcsr ) ) );
    }

    /// Computes y = alpha * A * x + beta * y
    HYPRE_Int Mult( HypreParVector& x, HypreParVector& y, double alpha = 1.0, double beta = 0.0 );
    /// Computes y = alpha * A * x + beta * y
    HYPRE_Int Mult( HYPRE_ParVector x, HYPRE_ParVector y, double alpha = 1.0, double beta = 0.0 );
    /// Computes y = alpha * A^t * x + beta * y
    HYPRE_Int MultTranspose( HypreParVector& x, HypreParVector& y, double alpha = 1.0, double beta = 0.0 );

    void Mult( double a, const HypreParVector& x, double b, HypreParVector& y ) const;
    void MultTranspose( double a, const HypreParVector& x, double b, HypreParVector& y ) const;

    // virtual void Mult(const Vector &x, Vector &y) const
    // { Mult(1.0, x, 0.0, y); }
    // virtual void MultTranspose(const Vector &x, Vector &y) const
    // { MultTranspose(1.0, x, 0.0, y); }

    /** The "Boolean" analog of y = alpha * A * x + beta * y, where elements in
        the sparsity pattern of the matrix are treated as "true". */
    // void BooleanMult(int alpha, int *x, int beta, int *y)
    // {
    //    internal::hypre_ParCSRMatrixBooleanMatvec(A_parcsr, alpha, x, beta, y);
    // }

    /** Multiply A on the left by a block-diagonal parallel matrix D. Return
        a new parallel matrix, D*A. If D has a different number of rows than A,
        D's row starts array needs to be given (as returned by the methods
        GetDofOffsets/GetTrueDofOffsets of ParFiniteElementSpace). The new matrix
        D*A uses copies of the row-, column-starts arrays, so "this" matrix and
        "row_starts" can be deleted.
        NOTE: this operation is local and does not require communication. */
    // HypreParMatrix* LeftDiagMult(const MOABSparseMatrix &D,
    //                              HYPRE_Int* row_starts = NULL) const;

    /// Scale the local row i by s(i).
    void ScaleRows( const Vector& s );
    /// Scale the local row i by 1./s(i)
    void InvScaleRows( const Vector& s );
    /// Scale all entries by s: A_scaled = s*A.
    void operator*=( double s );

    /// Remove values smaller in absolute value than some threshold
    void Threshold( double threshold = 0.0 );

    /// If a row contains only zeros, set its diagonal to 1.
    void EliminateZeroRows()
    {
        hypre_ParCSRMatrixFixZeroRows( A_parcsr );
    }

    /** Eliminate rows and columns from the matrix, and rows from the vector B.
        Modify B with the BC values in X. */
    void EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols, const HypreParVector& X, HypreParVector& B );

    /** Eliminate rows and columns from the matrix and store the eliminated
        elements in a new matrix Ae (returned), so that the modified matrix and
        Ae sum to the original matrix. */
    HypreParMatrix* EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols );

    HYPRE_Int GetValues( const HYPRE_Int nrows,
                         HYPRE_Int* ncols,
                         HYPRE_Int* rows,
                         HYPRE_Int* cols,
                         HYPRE_Complex* values );
    HYPRE_Int SetValues( const HYPRE_Int nrows,
                         HYPRE_Int* ncols,
                         const HYPRE_Int* rows,
                         const HYPRE_Int* cols,
                         const HYPRE_Complex* values );
    HYPRE_Int AddToValues( const HYPRE_Int nrows,
                           HYPRE_Int* ncols,
                           const HYPRE_Int* rows,
                           const HYPRE_Int* cols,
                           const HYPRE_Complex* values );

    HYPRE_Int FinalizeAssembly();

    HYPRE_Int verbosity( const HYPRE_Int level );

    /// Prints the locally owned rows in parallel
    void Print( const char* fname, HYPRE_Int offi = 0, HYPRE_Int offj = 0 );
    /// Reads the matrix from a file
    void Read( const char* fname );

    /// Calls hypre's destroy function
    virtual ~HypreParMatrix()
    {
        Destroy();
    }

    /// Returns the matrix A * B
    friend HypreParMatrix* ParMult( HypreParMatrix* A, HypreParMatrix* B );

    /// Returns the matrix P^t * A * P
    friend HypreParMatrix* RAP( HypreParMatrix* A, HypreParMatrix* P );
    /// Returns the matrix Rt^t * A * P
    friend HypreParMatrix* RAP( HypreParMatrix* Rt, HypreParMatrix* A, HypreParMatrix* P );

    /** Eliminate essential BC specified by 'ess_dof_list' from the solution X to
        the r.h.s. B. Here A is a matrix with eliminated BC, while Ae is such that
        (A+Ae) is the original (Neumann) matrix before elimination. */
    friend void EliminateBC( HypreParMatrix& A,
                             HypreParMatrix& Ae,
                             const std::vector< int >& ess_dof_list,
                             const HypreParVector& X,
                             HypreParVector& B );

    friend class HypreSolver;
};

}  // namespace moab

#endif  // MOAB_HAVE_MPI

#endif