MOAB: Mesh Oriented datABase  (version 5.2.1)
HypreParMatrix.hpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.org.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 /// This HYPRE library interface has been taken originally from MFEM and modified
00013 /// to suit the needs for the MOAB library.
00014 /// Modified by: Vijay Mahadevan
00015 
00016 #ifndef MOAB_HYPREPARMATRIX
00017 #define MOAB_HYPREPARMATRIX
00018 
00019 #include "moab/MOABConfig.h"
00020 #include "moab/Core.hpp"
00021 
00022 #ifdef MOAB_HAVE_EIGEN3
00023 #include <Eigen/Core>
00024 #include <Eigen/Sparse>
00025 #else
00026 #error Configure with Eigen3 enabled
00027 #endif
00028 
00029 #ifdef MOAB_HAVE_MPI
00030 #include "moab/ParallelComm.hpp"
00031 
00032 // hypre header files
00033 // #include "HYPRE.h"
00034 // #include "HYPRE_parcsr_ls.h"
00035 // #include "HYPRE_MatvecFunctions.h"
00036 
00037 #include "seq_mv.h"
00038 #include "HypreParVector.hpp"
00039 #include "_hypre_parcsr_mv.h"
00040 #include "_hypre_parcsr_ls.h"
00041 #include "interpreter.h"
00042 #include "HYPRE_MatvecFunctions.h"
00043 
00044 #include "temp_multivector.h"
00045 
00046 #ifdef HYPRE_COMPLEX
00047 #error "MOAB does not work with HYPRE's complex numbers support"
00048 #endif
00049 
00050 #include "hypre_parcsr.hpp"
00051 
00052 namespace moab
00053 {
00054 
00055 namespace internal
00056 {
00057 
00058     // Convert a HYPRE_Int to int
00059     inline int to_int( HYPRE_Int i )
00060     {
00061 #ifdef HYPRE_BIGINT
00062         MFEM_ASSERT( HYPRE_Int( int( i ) ) == i, "overflow converting HYPRE_Int to int" );
00063 #endif
00064         return int( i );
00065     }
00066 
00067 }  // namespace internal
00068 
00069 /// Wrapper for hypre's ParCSR matrix class
00070 class HypreParMatrix
00071 {
00072   private:
00073     /// The actual object
00074     HYPRE_IJMatrix A;
00075     hypre_ParCSRMatrix* A_parcsr;
00076 
00077     mv_InterfaceInterpreter* interpreter;
00078 
00079     // Does the object own the pointer A?
00080     char ParCSROwner;
00081 
00082     // Store the dimensions of the matrix
00083     int height, gnrows;
00084     int width, gncols;
00085 
00086     moab::ParallelComm* pcomm;
00087     char initialized;
00088 
00089     // Initialize with defaults. Does not initialize inherited members.
00090     void Init();
00091 
00092     // Delete all owned data. Does not perform re-initialization with defaults.
00093     void Destroy();
00094 
00095     friend class HypreSolver;
00096 
00097   public:
00098     typedef Eigen::VectorXd Vector;
00099     typedef Eigen::SparseMatrix< double, Eigen::RowMajor > MOABSparseMatrix;
00100     // typedef template Eigen::ArrayXX<T> Array2D<T>;
00101 
00102   public:
00103     /// An empty matrix to be used as a reference to an existing matrix
00104     HypreParMatrix( moab::ParallelComm* p_comm );
00105 
00106     /// Converts hypre's format to HypreParMatrix
00107     HypreParMatrix( HYPRE_IJMatrix a )
00108     {
00109         Init();
00110         A      = a;
00111         height = GetNumRows();
00112         width  = GetNumCols();
00113     }
00114 
00115     /** Creates block-diagonal square parallel matrix. Diagonal is given by diag
00116         which must be in CSR format (finalized). The new HypreParMatrix does not
00117         take ownership of any of the input arrays. */
00118     HypreParMatrix( moab::ParallelComm* p_comm, HYPRE_Int glob_size, HYPRE_Int* row_starts, HYPRE_Int nnz_pr_diag = 0 );
00119 
00120     /** Creates block-diagonal rectangular parallel matrix. Diagonal is given by
00121         diag which must be in CSR format (finalized). The new HypreParMatrix does
00122         not take ownership of any of the input arrays. */
00123     HypreParMatrix( moab::ParallelComm* p_comm, HYPRE_Int global_num_rows, HYPRE_Int global_num_cols,
00124                     HYPRE_Int* row_starts, HYPRE_Int* col_starts, HYPRE_Int nnz_pr_diag = 0, HYPRE_Int onz_pr_diag = 0,
00125                     HYPRE_Int nnz_pr_offdiag = 0 );
00126 
00127     /** Creates a general parallel matrix from a local CSR matrix on each
00128         processor described by the I, J and data arrays. The local matrix should
00129         be of size (local) nrows by (global) glob_ncols. The new parallel matrix
00130         contains copies of all input arrays (so they can be deleted). */
00131     // HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_Int glob_nrows,
00132     //                HYPRE_Int glob_ncols, int *I, HYPRE_Int *J,
00133     //                double *data, HYPRE_Int *rows, HYPRE_Int *cols);
00134 
00135     /// Make this HypreParMatrix a reference to 'master'
00136     void MakeRef( const HypreParMatrix& master );
00137 
00138     /// MPI communicator
00139     moab::ParallelComm* GetParallelCommunicator() const
00140     {
00141         return pcomm;
00142     }
00143 
00144     /// Typecasting to hypre's HYPRE_IJMatrix*
00145     operator HYPRE_IJMatrix()
00146     {
00147         return A;
00148     }
00149     /// Typecasting to hypre's HYPRE_ParCSRMatrix, a.k.a. void *
00150     operator HYPRE_ParCSRMatrix()
00151     {
00152         return (HYPRE_ParCSRMatrix)A_parcsr;
00153     }
00154 
00155     /// Changes the ownership of the matrix
00156     HYPRE_IJMatrix StealData();
00157     /// Changes the ownership of the matrix to A
00158     void StealData( HypreParMatrix& A );
00159 
00160     /** If the HypreParMatrix does not own the row-starts array, make a copy of
00161         it that the HypreParMatrix will own. If the col-starts array is the same
00162         as the row-starts array, col-starts is also replaced. */
00163     // void CopyRowStarts();
00164     /** If the HypreParMatrix does not own the col-starts array, make a copy of
00165         it that the HypreParMatrix will own. If the row-starts array is the same
00166         as the col-starts array, row-starts is also replaced. */
00167     // void CopyColStarts();
00168 
00169     /// Returns the global number of nonzeros
00170     inline HYPRE_Int NNZ()
00171     {
00172         return A_parcsr->num_nonzeros;
00173     }
00174     /// Returns the row partitioning
00175     inline HYPRE_Int* RowPart()
00176     {
00177         return A->row_partitioning;
00178     }
00179     /// Returns the column partitioning
00180     inline HYPRE_Int* ColPart()
00181     {
00182         return A->col_partitioning;
00183     }
00184     /// Returns the global number of rows
00185     inline HYPRE_Int M()
00186     {
00187         return A->global_num_rows;
00188     }
00189     /// Returns the global number of columns
00190     inline HYPRE_Int N()
00191     {
00192         return A->global_num_cols;
00193     }
00194     /// Returns the global number of rows
00195     inline HYPRE_Int FirstRow()
00196     {
00197         return A->global_first_row;
00198     }
00199     /// Returns the global number of columns
00200     inline HYPRE_Int FirstCol()
00201     {
00202         return A->global_first_col;
00203     }
00204 
00205     /// Get the local diagonal of the matrix.
00206     void GetDiag( Vector& diag ) const;
00207     /// Get the local diagonal block. NOTE: 'diag' will not own any data.
00208     void GetDiag( MOABSparseMatrix& diag ) const;
00209     /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
00210     void GetOffd( MOABSparseMatrix& offd, HYPRE_Int*& cmap ) const;
00211 
00212     /** Split the matrix into M x N equally sized blocks of parallel matrices.
00213         The size of 'blocks' must already be set to M x N. */
00214     // void GetBlocks(Eigen::Array<HypreParMatrix*,Eigen::Dynamic,Eigen::Dynamic> &blocks,
00215     //                bool interleaved_rows = false,
00216     //                bool interleaved_cols = false) const;
00217 
00218     /// Returns the transpose of *this
00219     // HypreParMatrix * Transpose();
00220 
00221     /** Destroy and resize block-diagonal square parallel matrix. Diagonal is given by diag
00222         which must be in CSR format (finalized). The new HypreParMatrix does not
00223         take ownership of any of the input arrays. */
00224     void resize( HYPRE_Int glob_size, HYPRE_Int* row_starts, HYPRE_Int nnz_pr_diag = 0, HYPRE_Int nnz_pr_offdiag = 0 );
00225 
00226     /** Destroy and resize block-diagonal rectangular parallel matrix. Diagonal is given by
00227         diag which must be in CSR format (finalized). The new HypreParMatrix does
00228         not take ownership of any of the input arrays. */
00229     void resize( HYPRE_Int global_num_rows, HYPRE_Int global_num_cols, HYPRE_Int* row_starts, HYPRE_Int* col_starts,
00230                  HYPRE_Int* nnz_pr_diag = NULL, HYPRE_Int* onz_pr_diag = NULL, HYPRE_Int nnz_pr_offdiag = 0 );
00231 
00232     /// Returns the number of rows in the diagonal block of the ParCSRMatrix
00233     int GetNumRows() const
00234     {
00235         return internal::to_int( hypre_CSRMatrixNumRows( hypre_ParCSRMatrixDiag( A_parcsr ) ) );
00236     }
00237 
00238     /// Returns the number of columns in the diagonal block of the ParCSRMatrix
00239     int GetNumCols() const
00240     {
00241         return internal::to_int( hypre_CSRMatrixNumCols( hypre_ParCSRMatrixDiag( A_parcsr ) ) );
00242     }
00243 
00244     /// Computes y = alpha * A * x + beta * y
00245     HYPRE_Int Mult( HypreParVector& x, HypreParVector& y, double alpha = 1.0, double beta = 0.0 );
00246     /// Computes y = alpha * A * x + beta * y
00247     HYPRE_Int Mult( HYPRE_ParVector x, HYPRE_ParVector y, double alpha = 1.0, double beta = 0.0 );
00248     /// Computes y = alpha * A^t * x + beta * y
00249     HYPRE_Int MultTranspose( HypreParVector& x, HypreParVector& y, double alpha = 1.0, double beta = 0.0 );
00250 
00251     void Mult( double a, const HypreParVector& x, double b, HypreParVector& y ) const;
00252     void MultTranspose( double a, const HypreParVector& x, double b, HypreParVector& y ) const;
00253 
00254     // virtual void Mult(const Vector &x, Vector &y) const
00255     // { Mult(1.0, x, 0.0, y); }
00256     // virtual void MultTranspose(const Vector &x, Vector &y) const
00257     // { MultTranspose(1.0, x, 0.0, y); }
00258 
00259     /** The "Boolean" analog of y = alpha * A * x + beta * y, where elements in
00260         the sparsity pattern of the matrix are treated as "true". */
00261     // void BooleanMult(int alpha, int *x, int beta, int *y)
00262     // {
00263     //    internal::hypre_ParCSRMatrixBooleanMatvec(A_parcsr, alpha, x, beta, y);
00264     // }
00265 
00266     /** Multiply A on the left by a block-diagonal parallel matrix D. Return
00267         a new parallel matrix, D*A. If D has a different number of rows than A,
00268         D's row starts array needs to be given (as returned by the methods
00269         GetDofOffsets/GetTrueDofOffsets of ParFiniteElementSpace). The new matrix
00270         D*A uses copies of the row-, column-starts arrays, so "this" matrix and
00271         "row_starts" can be deleted.
00272         NOTE: this operation is local and does not require communication. */
00273     // HypreParMatrix* LeftDiagMult(const MOABSparseMatrix &D,
00274     //                              HYPRE_Int* row_starts = NULL) const;
00275 
00276     /// Scale the local row i by s(i).
00277     void ScaleRows( const Vector& s );
00278     /// Scale the local row i by 1./s(i)
00279     void InvScaleRows( const Vector& s );
00280     /// Scale all entries by s: A_scaled = s*A.
00281     void operator*=( double s );
00282 
00283     /// Remove values smaller in absolute value than some threshold
00284     void Threshold( double threshold = 0.0 );
00285 
00286     /// If a row contains only zeros, set its diagonal to 1.
00287     void EliminateZeroRows()
00288     {
00289         hypre_ParCSRMatrixFixZeroRows( A_parcsr );
00290     }
00291 
00292     /** Eliminate rows and columns from the matrix, and rows from the vector B.
00293         Modify B with the BC values in X. */
00294     void EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols, const HypreParVector& X, HypreParVector& B );
00295 
00296     /** Eliminate rows and columns from the matrix and store the eliminated
00297         elements in a new matrix Ae (returned), so that the modified matrix and
00298         Ae sum to the original matrix. */
00299     HypreParMatrix* EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols );
00300 
00301     HYPRE_Int GetValues( const HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_Int* rows, HYPRE_Int* cols,
00302                          HYPRE_Complex* values );
00303     HYPRE_Int SetValues( const HYPRE_Int nrows, HYPRE_Int* ncols, const HYPRE_Int* rows, const HYPRE_Int* cols,
00304                          const HYPRE_Complex* values );
00305     HYPRE_Int AddToValues( const HYPRE_Int nrows, HYPRE_Int* ncols, const HYPRE_Int* rows, const HYPRE_Int* cols,
00306                            const HYPRE_Complex* values );
00307 
00308     HYPRE_Int FinalizeAssembly();
00309 
00310     HYPRE_Int verbosity( const HYPRE_Int level );
00311 
00312     /// Prints the locally owned rows in parallel
00313     void Print( const char* fname, HYPRE_Int offi = 0, HYPRE_Int offj = 0 );
00314     /// Reads the matrix from a file
00315     void Read( const char* fname );
00316 
00317     /// Calls hypre's destroy function
00318     virtual ~HypreParMatrix()
00319     {
00320         Destroy();
00321     }
00322 
00323     /// Returns the matrix A * B
00324     friend HypreParMatrix* ParMult( HypreParMatrix* A, HypreParMatrix* B );
00325 
00326     /// Returns the matrix P^t * A * P
00327     friend HypreParMatrix* RAP( HypreParMatrix* A, HypreParMatrix* P );
00328     /// Returns the matrix Rt^t * A * P
00329     friend HypreParMatrix* RAP( HypreParMatrix* Rt, HypreParMatrix* A, HypreParMatrix* P );
00330 
00331     /** Eliminate essential BC specified by 'ess_dof_list' from the solution X to
00332         the r.h.s. B. Here A is a matrix with eliminated BC, while Ae is such that
00333         (A+Ae) is the original (Neumann) matrix before elimination. */
00334     friend void EliminateBC( HypreParMatrix& A, HypreParMatrix& Ae, const std::vector< int >& ess_dof_list,
00335                              const HypreParVector& X, HypreParVector& B );
00336 
00337     friend class HypreSolver;
00338 };
00339 
00340 }  // namespace moab
00341 
00342 #endif  // MOAB_HAVE_MPI
00343 
00344 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines