MOAB: Mesh Oriented datABase  (version 5.4.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,
00124                     HYPRE_Int global_num_rows,
00125                     HYPRE_Int global_num_cols,
00126                     HYPRE_Int* row_starts,
00127                     HYPRE_Int* col_starts,
00128                     HYPRE_Int nnz_pr_diag    = 0,
00129                     HYPRE_Int onz_pr_diag    = 0,
00130                     HYPRE_Int nnz_pr_offdiag = 0 );
00131 
00132     /** Creates a general parallel matrix from a local CSR matrix on each
00133         processor described by the I, J and data arrays. The local matrix should
00134         be of size (local) nrows by (global) glob_ncols. The new parallel matrix
00135         contains copies of all input arrays (so they can be deleted). */
00136     // HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_Int glob_nrows,
00137     //                HYPRE_Int glob_ncols, int *I, HYPRE_Int *J,
00138     //                double *data, HYPRE_Int *rows, HYPRE_Int *cols);
00139 
00140     /// Make this HypreParMatrix a reference to 'master'
00141     void MakeRef( const HypreParMatrix& master );
00142 
00143     /// MPI communicator
00144     moab::ParallelComm* GetParallelCommunicator() const
00145     {
00146         return pcomm;
00147     }
00148 
00149     /// Typecasting to hypre's HYPRE_IJMatrix*
00150     operator HYPRE_IJMatrix()
00151     {
00152         return A;
00153     }
00154     /// Typecasting to hypre's HYPRE_ParCSRMatrix, a.k.a. void *
00155     operator HYPRE_ParCSRMatrix()
00156     {
00157         return (HYPRE_ParCSRMatrix)A_parcsr;
00158     }
00159 
00160     /// Changes the ownership of the matrix
00161     HYPRE_IJMatrix StealData();
00162     /// Changes the ownership of the matrix to A
00163     void StealData( HypreParMatrix& A );
00164 
00165     /** If the HypreParMatrix does not own the row-starts array, make a copy of
00166         it that the HypreParMatrix will own. If the col-starts array is the same
00167         as the row-starts array, col-starts is also replaced. */
00168     // void CopyRowStarts();
00169     /** If the HypreParMatrix does not own the col-starts array, make a copy of
00170         it that the HypreParMatrix will own. If the row-starts array is the same
00171         as the col-starts array, row-starts is also replaced. */
00172     // void CopyColStarts();
00173 
00174     /// Returns the global number of nonzeros
00175     inline HYPRE_Int NNZ()
00176     {
00177         return A_parcsr->num_nonzeros;
00178     }
00179     /// Returns the row partitioning
00180     inline HYPRE_Int* RowPart()
00181     {
00182         return A->row_partitioning;
00183     }
00184     /// Returns the column partitioning
00185     inline HYPRE_Int* ColPart()
00186     {
00187         return A->col_partitioning;
00188     }
00189     /// Returns the global number of rows
00190     inline HYPRE_Int M()
00191     {
00192         return A->global_num_rows;
00193     }
00194     /// Returns the global number of columns
00195     inline HYPRE_Int N()
00196     {
00197         return A->global_num_cols;
00198     }
00199     /// Returns the global number of rows
00200     inline HYPRE_Int FirstRow()
00201     {
00202         return A->global_first_row;
00203     }
00204     /// Returns the global number of columns
00205     inline HYPRE_Int FirstCol()
00206     {
00207         return A->global_first_col;
00208     }
00209 
00210     /// Get the local diagonal of the matrix.
00211     void GetDiag( Vector& diag ) const;
00212     /// Get the local diagonal block. NOTE: 'diag' will not own any data.
00213     void GetDiag( MOABSparseMatrix& diag ) const;
00214     /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
00215     void GetOffd( MOABSparseMatrix& offd, HYPRE_Int*& cmap ) const;
00216 
00217     /** Split the matrix into M x N equally sized blocks of parallel matrices.
00218         The size of 'blocks' must already be set to M x N. */
00219     // void GetBlocks(Eigen::Array<HypreParMatrix*,Eigen::Dynamic,Eigen::Dynamic> &blocks,
00220     //                bool interleaved_rows = false,
00221     //                bool interleaved_cols = false) const;
00222 
00223     /// Returns the transpose of *this
00224     // HypreParMatrix * Transpose();
00225 
00226     /** Destroy and resize block-diagonal square parallel matrix. Diagonal is given by diag
00227         which must be in CSR format (finalized). The new HypreParMatrix does not
00228         take ownership of any of the input arrays. */
00229     void resize( HYPRE_Int glob_size, HYPRE_Int* row_starts, HYPRE_Int nnz_pr_diag = 0, HYPRE_Int nnz_pr_offdiag = 0 );
00230 
00231     /** Destroy and resize block-diagonal rectangular parallel matrix. Diagonal is given by
00232         diag which must be in CSR format (finalized). The new HypreParMatrix does
00233         not take ownership of any of the input arrays. */
00234     void resize( HYPRE_Int global_num_rows,
00235                  HYPRE_Int global_num_cols,
00236                  HYPRE_Int* row_starts,
00237                  HYPRE_Int* col_starts,
00238                  HYPRE_Int* nnz_pr_diag   = NULL,
00239                  HYPRE_Int* onz_pr_diag   = NULL,
00240                  HYPRE_Int nnz_pr_offdiag = 0 );
00241 
00242     /// Returns the number of rows in the diagonal block of the ParCSRMatrix
00243     int GetNumRows() const
00244     {
00245         return internal::to_int( hypre_CSRMatrixNumRows( hypre_ParCSRMatrixDiag( A_parcsr ) ) );
00246     }
00247 
00248     /// Returns the number of columns in the diagonal block of the ParCSRMatrix
00249     int GetNumCols() const
00250     {
00251         return internal::to_int( hypre_CSRMatrixNumCols( hypre_ParCSRMatrixDiag( A_parcsr ) ) );
00252     }
00253 
00254     /// Computes y = alpha * A * x + beta * y
00255     HYPRE_Int Mult( HypreParVector& x, HypreParVector& y, double alpha = 1.0, double beta = 0.0 );
00256     /// Computes y = alpha * A * x + beta * y
00257     HYPRE_Int Mult( HYPRE_ParVector x, HYPRE_ParVector y, double alpha = 1.0, double beta = 0.0 );
00258     /// Computes y = alpha * A^t * x + beta * y
00259     HYPRE_Int MultTranspose( HypreParVector& x, HypreParVector& y, double alpha = 1.0, double beta = 0.0 );
00260 
00261     void Mult( double a, const HypreParVector& x, double b, HypreParVector& y ) const;
00262     void MultTranspose( double a, const HypreParVector& x, double b, HypreParVector& y ) const;
00263 
00264     // virtual void Mult(const Vector &x, Vector &y) const
00265     // { Mult(1.0, x, 0.0, y); }
00266     // virtual void MultTranspose(const Vector &x, Vector &y) const
00267     // { MultTranspose(1.0, x, 0.0, y); }
00268 
00269     /** The "Boolean" analog of y = alpha * A * x + beta * y, where elements in
00270         the sparsity pattern of the matrix are treated as "true". */
00271     // void BooleanMult(int alpha, int *x, int beta, int *y)
00272     // {
00273     //    internal::hypre_ParCSRMatrixBooleanMatvec(A_parcsr, alpha, x, beta, y);
00274     // }
00275 
00276     /** Multiply A on the left by a block-diagonal parallel matrix D. Return
00277         a new parallel matrix, D*A. If D has a different number of rows than A,
00278         D's row starts array needs to be given (as returned by the methods
00279         GetDofOffsets/GetTrueDofOffsets of ParFiniteElementSpace). The new matrix
00280         D*A uses copies of the row-, column-starts arrays, so "this" matrix and
00281         "row_starts" can be deleted.
00282         NOTE: this operation is local and does not require communication. */
00283     // HypreParMatrix* LeftDiagMult(const MOABSparseMatrix &D,
00284     //                              HYPRE_Int* row_starts = NULL) const;
00285 
00286     /// Scale the local row i by s(i).
00287     void ScaleRows( const Vector& s );
00288     /// Scale the local row i by 1./s(i)
00289     void InvScaleRows( const Vector& s );
00290     /// Scale all entries by s: A_scaled = s*A.
00291     void operator*=( double s );
00292 
00293     /// Remove values smaller in absolute value than some threshold
00294     void Threshold( double threshold = 0.0 );
00295 
00296     /// If a row contains only zeros, set its diagonal to 1.
00297     void EliminateZeroRows()
00298     {
00299         hypre_ParCSRMatrixFixZeroRows( A_parcsr );
00300     }
00301 
00302     /** Eliminate rows and columns from the matrix, and rows from the vector B.
00303         Modify B with the BC values in X. */
00304     void EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols, const HypreParVector& X, HypreParVector& B );
00305 
00306     /** Eliminate rows and columns from the matrix and store the eliminated
00307         elements in a new matrix Ae (returned), so that the modified matrix and
00308         Ae sum to the original matrix. */
00309     HypreParMatrix* EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols );
00310 
00311     HYPRE_Int GetValues( const HYPRE_Int nrows,
00312                          HYPRE_Int* ncols,
00313                          HYPRE_Int* rows,
00314                          HYPRE_Int* cols,
00315                          HYPRE_Complex* values );
00316     HYPRE_Int SetValues( const HYPRE_Int nrows,
00317                          HYPRE_Int* ncols,
00318                          const HYPRE_Int* rows,
00319                          const HYPRE_Int* cols,
00320                          const HYPRE_Complex* values );
00321     HYPRE_Int AddToValues( const HYPRE_Int nrows,
00322                            HYPRE_Int* ncols,
00323                            const HYPRE_Int* rows,
00324                            const HYPRE_Int* cols,
00325                            const HYPRE_Complex* values );
00326 
00327     HYPRE_Int FinalizeAssembly();
00328 
00329     HYPRE_Int verbosity( const HYPRE_Int level );
00330 
00331     /// Prints the locally owned rows in parallel
00332     void Print( const char* fname, HYPRE_Int offi = 0, HYPRE_Int offj = 0 );
00333     /// Reads the matrix from a file
00334     void Read( const char* fname );
00335 
00336     /// Calls hypre's destroy function
00337     virtual ~HypreParMatrix()
00338     {
00339         Destroy();
00340     }
00341 
00342     /// Returns the matrix A * B
00343     friend HypreParMatrix* ParMult( HypreParMatrix* A, HypreParMatrix* B );
00344 
00345     /// Returns the matrix P^t * A * P
00346     friend HypreParMatrix* RAP( HypreParMatrix* A, HypreParMatrix* P );
00347     /// Returns the matrix Rt^t * A * P
00348     friend HypreParMatrix* RAP( HypreParMatrix* Rt, HypreParMatrix* A, HypreParMatrix* P );
00349 
00350     /** Eliminate essential BC specified by 'ess_dof_list' from the solution X to
00351         the r.h.s. B. Here A is a matrix with eliminated BC, while Ae is such that
00352         (A+Ae) is the original (Neumann) matrix before elimination. */
00353     friend void EliminateBC( HypreParMatrix& A,
00354                              HypreParMatrix& Ae,
00355                              const std::vector< int >& ess_dof_list,
00356                              const HypreParVector& X,
00357                              HypreParVector& B );
00358 
00359     friend class HypreSolver;
00360 };
00361 
00362 }  // namespace moab
00363 
00364 #endif  // MOAB_HAVE_MPI
00365 
00366 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines