HypreParMatrix.cpp

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.
 
#include "HypreParVector.hpp"
#include "HypreParMatrix.hpp"
 
#ifdef MOAB_HAVE_MPI
 
#include "hypre_parcsr.hpp"
 
#include <fstream>
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cassert>
 
using namespace std;
 
#define MOAB_DEBUG
 
namespace moab
{
#define moab_hypre_assert( a, b ) \
 { \
 }
#define moab_hypre_assert_t( a, b ) \
 { \
 if( !a ) \
 { \
 std::cout << "HYPRE Error: " << b << std::endl; \
 exit( -1 ); \
 } \
 }
 
template < typename TargetT, typename SourceT >
static TargetT* DuplicateAs( const SourceT* array, int size, bool cplusplus = true )
{
 TargetT* target_array = cplusplus ? new TargetT[size]
 /* */
 : hypre_TAlloc( TargetT, size );
 
 for( int i = 0; i < size; i++ )
 {
 target_array[i] = array[i];
 }
 
 return target_array;
}
 
void HypreParMatrix::Init()
{
 A = NULL;
 A_parcsr = NULL;
 ParCSROwner = 1;
}
 
HypreParMatrix::HypreParMatrix( moab::ParallelComm* p_comm ) : pcomm( p_comm )
{
 Init();
 height = width = 0;
}
 
// Square block-diagonal constructor
HypreParMatrix::HypreParMatrix( moab::ParallelComm* p_comm,
 HYPRE_Int glob_size,
 HYPRE_Int* row_starts,
 HYPRE_Int nnz_pr )
 : pcomm( p_comm )
{
 Init();
 resize( glob_size, row_starts, nnz_pr );
}
 
void HypreParMatrix::resize( HYPRE_Int glob_size,
 HYPRE_Int* row_starts,
 HYPRE_Int nnz_pr_diag,
 HYPRE_Int nnz_pr_offdiag )
{
 /* Create the matrix.
 Note that this is a square matrix, so we indicate the row partition
 size twice (since number of rows = number of cols) */
 HYPRE_IJMatrixCreate( pcomm->comm(), row_starts[0], row_starts[1], row_starts[0], row_starts[1], &A );
 /* Choose a parallel csr format storage (see the User's Manual) */
 HYPRE_IJMatrixSetObjectType( A, HYPRE_PARCSR );
 
 int locsize = row_starts[1] - row_starts[0] + 1;
 int num_procs, rank;
 MPI_Comm_size( pcomm->comm(), &num_procs );
 MPI_Comm_rank( pcomm->comm(), &rank );
 
 if( nnz_pr_diag != 0 || nnz_pr_offdiag != 0 )
 {
 if( num_procs > 1 )
 {
 std::vector< HYPRE_Int > m_nnz_pr_diag( locsize, nnz_pr_diag );
 std::vector< HYPRE_Int > m_onz_pr_diag( locsize, nnz_pr_offdiag );
 HYPRE_IJMatrixSetDiagOffdSizes( A, &m_nnz_pr_diag[0], &m_onz_pr_diag[0] );
 }
 else
 {
 std::vector< HYPRE_Int > m_nnz_pr_diag( locsize, nnz_pr_diag );
 HYPRE_IJMatrixSetRowSizes( A, &m_nnz_pr_diag[0] );
 }
 }
 
 /* Initialize before setting coefficients */
 HYPRE_IJMatrixInitialize( A );
 /* Get the parcsr matrix object to use */
 HYPRE_IJMatrixGetObject( A, (void**)&A_parcsr );
 
 /* define an interpreter for the ParCSR interface */
 HYPRE_MatvecFunctions matvec_fn;
 interpreter = hypre_CTAlloc( mv_InterfaceInterpreter, 1 );
 HYPRE_ParCSRSetupInterpreter( interpreter );
 HYPRE_ParCSRSetupMatvec( &matvec_fn );
 gnrows = glob_size;
 gncols = glob_size;
 height = GetNumRows();
 width = GetNumCols();
}
 
void HypreParMatrix::resize( HYPRE_Int global_num_rows,
 HYPRE_Int global_num_cols,
 HYPRE_Int* row_starts,
 HYPRE_Int* col_starts,
 HYPRE_Int* nnz_pr_diag,
 HYPRE_Int* onz_pr_diag,
 HYPRE_Int nnz_pr_offdiag )
{
 /* Create the matrix.
 Note that this is a square matrix, so we indicate the row partition
 size twice (since number of rows = number of cols) */
 HYPRE_IJMatrixCreate( pcomm->comm(), row_starts[0], row_starts[1], col_starts[0], col_starts[1], &A );
 /* Choose a parallel csr format storage (see the User's Manual) */
 HYPRE_IJMatrixSetObjectType( A, HYPRE_PARCSR );
 
 int num_procs;
 MPI_Comm_size( pcomm->comm(), &num_procs );
 
 /* Set the matrix pre-allocation data */
 if( num_procs > 1 )
 {
 if( nnz_pr_diag == NULL && onz_pr_diag == NULL )
 {
 std::cout << "Parameter nnz_pr_diag and onz_pr_diag cannot be NULL.\n";
 moab_hypre_assert_t( nnz_pr_diag, true );
 }
 
 HYPRE_IJMatrixSetDiagOffdSizes( A, nnz_pr_diag, onz_pr_diag );
 }
 else
 {
 HYPRE_IJMatrixSetRowSizes( A, nnz_pr_diag );
 }
 
 if( nnz_pr_offdiag )
 {
 HYPRE_IJMatrixSetMaxOffProcElmts( A, nnz_pr_offdiag );
 }
 
 /* Initialize before setting coefficients */
 HYPRE_IJMatrixInitialize( A );
 /* Get the parcsr matrix object to use */
 HYPRE_IJMatrixGetObject( A, (void**)&A_parcsr );
 
 /* define an interpreter for the ParCSR interface */
 HYPRE_MatvecFunctions matvec_fn;
 interpreter = hypre_CTAlloc( mv_InterfaceInterpreter, 1 );
 HYPRE_ParCSRSetupInterpreter( interpreter );
 HYPRE_ParCSRSetupMatvec( &matvec_fn );
 gnrows = global_num_rows;
 gncols = global_num_cols;
 height = GetNumRows();
 width = GetNumCols();
}
 
// Rectangular block-diagonal constructor
HypreParMatrix::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,
 HYPRE_Int onz_pr_diag,
 HYPRE_Int nnz_pr_offdiag )
 : pcomm( p_comm )
{
 Init();
 std::vector< HYPRE_Int > m_nnz_pr_diag( row_starts[1] - row_starts[0], nnz_pr_diag );
 std::vector< HYPRE_Int > m_onz_pr_diag( row_starts[1] - row_starts[0], onz_pr_diag );
 resize( global_num_rows, global_num_cols, row_starts, col_starts, &m_nnz_pr_diag[0], &m_onz_pr_diag[0],
 nnz_pr_offdiag );
}
 
void HypreParMatrix::MakeRef( const HypreParMatrix& master )
{
 Destroy();
 Init();
 A = master.A;
 A_parcsr = master.A_parcsr;
 ParCSROwner = 0;
 height = master.GetNumRows();
 width = master.GetNumCols();
}
 
HYPRE_IJMatrix HypreParMatrix::StealData()
{
 // Only safe when (diagOwner == -1 && offdOwner == -1 && colMapOwner == -1)
 // Otherwise, there may be memory leaks or hypre may destroy arrays allocated
 // with operator new.
 moab_hypre_assert( diagOwner == -1 && offdOwner == -1 && colMapOwner == -1, "" );
 moab_hypre_assert( ParCSROwner, "" );
 HYPRE_IJMatrix R = A;
 A = NULL;
 ParCSROwner = 0;
 Destroy();
 Init();
 return R;
}
 
void HypreParMatrix::StealData( HypreParMatrix& X )
{
 // Only safe when (diagOwner == -1 && offdOwner == -1 && colMapOwner == -1)
 // Otherwise, there may be memory leaks or hypre may destroy arrays allocated
 // with operator new.
 moab_hypre_assert( diagOwner == -1 && offdOwner == -1 && colMapOwner == -1, "" );
 moab_hypre_assert( ParCSROwner, "" );
 moab_hypre_assert( X.diagOwner == -1 && X.offdOwner == -1 && X.colMapOwner == -1, "" );
 moab_hypre_assert( X.ParCSROwner, "" );
 A = X.A;
 A_parcsr = X.A_parcsr;
 ParCSROwner = 1;
 X.A = NULL;
 X.A_parcsr = NULL;
 X.ParCSROwner = 0;
 X.Destroy();
 X.Init();
}
 
void HypreParMatrix::GetDiag( Eigen::VectorXd& diag ) const
{
 int size = this->height;
 diag.resize( size );
 
 for( int j = 0; j < size; j++ )
 {
 diag( j ) = A_parcsr->diag->data[A_parcsr->diag->i[j]];
 moab_hypre_assert( A_parcsr->diag->j[A_parcsr->diag->i[j]] == j,
 "the first entry in each row must be the diagonal one" );
 }
}
 
static void MakeWrapper( const hypre_CSRMatrix* mat, Eigen::SparseMatrix< double, Eigen::RowMajor >& wrapper )
{
 HYPRE_Int nr = hypre_CSRMatrixNumRows( mat );
 HYPRE_Int nc = hypre_CSRMatrixNumCols( mat );
 Eigen::SparseMatrix< double, Eigen::RowMajor > tmp( nr, nc );
 typedef Eigen::Triplet< double > T;
 HYPRE_Int nnz = hypre_CSRMatrixNumNonzeros( mat );
 int *rindx, *cindx;
 double* dindx;
#ifndef HYPRE_BIGINT
 rindx = hypre_CSRMatrixI( mat );
 cindx = hypre_CSRMatrixJ( mat );
 dindx = hypre_CSRMatrixData( mat );
#else
 rindx = DuplicateAs< int >( hypre_CSRMatrixI( mat ), nr + 1 );
 cindx = DuplicateAs< int >( hypre_CSRMatrixJ( mat ), nnz );
 dindx = hypre_CSRMatrixData( mat );
#endif
 std::vector< T > tripletList( nnz );
 
 for( int i = 0; i < nnz; ++i )
 {
 tripletList.push_back( T( rindx[i], cindx[i], dindx[i] ) );
 }
 
 tmp.setFromTriplets( tripletList.begin(), tripletList.end() );
 wrapper.swap( tmp );
}
 
void HypreParMatrix::GetDiag( Eigen::SparseMatrix< double, Eigen::RowMajor >& diag ) const
{
 MakeWrapper( A_parcsr->diag, diag );
}
 
void HypreParMatrix::GetOffd( Eigen::SparseMatrix< double, Eigen::RowMajor >& offd, HYPRE_Int*& cmap ) const
{
 MakeWrapper( A_parcsr->offd, offd );
 cmap = A_parcsr->col_map_offd;
}
 
// void HypreParMatrix::GetBlocks(Eigen::Array<HypreParMatrix*,Eigen::Dynamic,Eigen::Dynamic>
// &blocks,
// bool interleaved_rows,
// bool interleaved_cols) const
// {
// int nr = blocks.rows();
// int nc = blocks.cols();
 
// hypre_ParCSRMatrix **hypre_blocks = new hypre_ParCSRMatrix*[nr * nc];
// internal::hypre_ParCSRMatrixSplit(A, nr, nc, hypre_blocks,
// interleaved_rows, interleaved_cols);
 
// for (int i = 0; i < nr; i++)
// {
// for (int j = 0; j < nc; j++)
// {
// blocks[i][j] = new HypreParMatrix(hypre_blocks[i*nc + j]);
// }
// }
 
// delete [] hypre_blocks;
// }
 
// HypreParMatrix * HypreParMatrix::Transpose()
// {
// hypre_ParCSRMatrix * At;
// hypre_ParCSRMatrixTranspose(A_parcsr, &At, 1);
// hypre_ParCSRMatrixSetNumNonzeros(At);
 
// hypre_MatvecCommPkgCreate(At);
 
// return new HypreParMatrix(At);
// }
 
HYPRE_Int HypreParMatrix::Mult( HypreParVector& x, HypreParVector& y, double a, double b )
{
 return hypre_ParCSRMatrixMatvec( a, A_parcsr, x.x_par, b, y.x_par );
}
 
void HypreParMatrix::Mult( double a, const HypreParVector& x, double b, HypreParVector& y ) const
{
 hypre_ParCSRMatrixMatvec( a, A_parcsr, x.x_par, b, y.x_par );
}
 
void HypreParMatrix::MultTranspose( double a, const HypreParVector& x, double b, HypreParVector& y ) const
{
 hypre_ParCSRMatrixMatvecT( a, A_parcsr, y.x_par, b, x.x_par );
}
 
HYPRE_Int HypreParMatrix::Mult( HYPRE_ParVector x, HYPRE_ParVector y, double a, double b )
{
 return hypre_ParCSRMatrixMatvec( a, A_parcsr, (hypre_ParVector*)x, b, (hypre_ParVector*)y );
}
 
HYPRE_Int HypreParMatrix::MultTranspose( HypreParVector& x, HypreParVector& y, double a, double b )
{
 return hypre_ParCSRMatrixMatvecT( a, A_parcsr, x.x_par, b, y.x_par );
}
 
// HypreParMatrix* HypreParMatrix::LeftDiagMult(const Eigen::SparseMatrix<double, Eigen::RowMajor>
// &D,
// HYPRE_Int* row_starts) const
// {
// const bool assumed_partition = HYPRE_AssumedPartitionCheck();
// const bool same_rows = (D.rows() == hypre_CSRMatrixNumRows(A_parcsr->diag));
 
// int part_size;
// if (assumed_partition)
// {
// part_size = 2;
// }
// else
// {
// MPI_Comm_size(GetComm(), &part_size);
// part_size++;
// }
 
// HYPRE_Int global_num_rows;
// if (same_rows)
// {
// row_starts = hypre_ParCSRMatrixRowStarts(A_parcsr);
// global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A_parcsr);
// }
// else
// {
// // MFEM_VERIFY(row_starts != NULL, "the number of rows in D and A is not "
// // "the same; row_starts must be given (not NULL)");
// // Here, when assumed_partition is true we use row_starts[2], so
// // row_starts must come from the GetDofOffsets/GetTrueDofOffsets methods
// // of ParFiniteElementSpace (HYPRE's partitions have only 2 entries).
// global_num_rows =
// assumed_partition ? row_starts[2] : row_starts[part_size-1];
// }
 
// HYPRE_Int *col_starts = hypre_ParCSRMatrixColStarts(A_parcsr);
// HYPRE_Int *col_map_offd;
 
// // get the diag and offd blocks as SparseMatrix wrappers
// Eigen::SparseMatrix<double, Eigen::RowMajor> A_diag, A_offd;
// GetDiag(A_diag);
// GetOffd(A_offd, col_map_offd);
 
// // multiply the blocks with D and create a new HypreParMatrix
// Eigen::SparseMatrix<double, Eigen::RowMajor> DA_diag = (D * A_diag);
// Eigen::SparseMatrix<double, Eigen::RowMajor> DA_offd = (D * A_offd);
 
// HypreParMatrix* DA =
// new HypreParMatrix(GetComm(),
// global_num_rows, hypre_ParCSRMatrixGlobalNumCols(A),
// DuplicateAs<HYPRE_Int>(row_starts, part_size, false),
// DuplicateAs<HYPRE_Int>(col_starts, part_size, false),
// &DA_diag, &DA_offd,
// DuplicateAs<HYPRE_Int>(col_map_offd, A_offd.cols()));
 
// // When HYPRE_BIGINT is defined, we want DA_{diag,offd} to delete their I and
// // J arrays but not their data arrays; when HYPRE_BIGINT is not defined, we
// // don't want DA_{diag,offd} to delete anything.
// // #ifndef HYPRE_BIGINT
// // DA_diag.LoseData();
// // DA_offd.LoseData();
// // #else
// // DA_diag.SetDataOwner(false);
// // DA_offd.SetDataOwner(false);
// // #endif
 
// // delete DA_diag;
// // delete DA_offd;
 
// hypre_ParCSRMatrixSetRowStartsOwner(DA->A, 1);
// hypre_ParCSRMatrixSetColStartsOwner(DA->A, 1);
 
// DA->diagOwner = DA->offdOwner = 3;
// DA->colMapOwner = 1;
 
// return DA;
// }
 
void HypreParMatrix::ScaleRows( const Eigen::VectorXd& diag )
{
 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
 {
 MB_SET_ERR_RET( "Row does not match" );
 }
 
 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != diag.size() )
 {
 MB_SET_ERR_RET( "Note the Eigen::VectorXd diag is not of compatible dimensions with A\n" );
 }
 
 int size = this->height;
 double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
 HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
 double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
 HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
 double val;
 HYPRE_Int jj;
 
 for( int i( 0 ); i < size; ++i )
 {
 val = diag[i];
 
 for( jj = Adiag_i[i]; jj < Adiag_i[i + 1]; ++jj )
 {
 Adiag_data[jj] *= val;
 }
 
 for( jj = Aoffd_i[i]; jj < Aoffd_i[i + 1]; ++jj )
 {
 Aoffd_data[jj] *= val;
 }
 }
}
 
void HypreParMatrix::InvScaleRows( const Eigen::VectorXd& diag )
{
 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
 {
 MB_SET_ERR_RET( "Row does not match" );
 }
 
 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != diag.size() )
 {
 MB_SET_ERR_RET( "Note the Eigen::VectorXd diag is not of compatible dimensions with A_parcsr\n" );
 }
 
 int size = this->height;
 double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
 HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
 double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
 HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
 double val;
 HYPRE_Int jj;
 
 for( int i( 0 ); i < size; ++i )
 {
#ifdef MOAB_DEBUG
 
 if( 0.0 == diag( i ) )
 {
 MB_SET_ERR_RET( "HypreParMatrix::InvDiagScale : Division by 0" );
 }
 
#endif
 val = 1. / diag( i );
 
 for( jj = Adiag_i[i]; jj < Adiag_i[i + 1]; ++jj )
 {
 Adiag_data[jj] *= val;
 }
 
 for( jj = Aoffd_i[i]; jj < Aoffd_i[i + 1]; ++jj )
 {
 Aoffd_data[jj] *= val;
 }
 }
}
 
void HypreParMatrix::operator*=( double s )
{
 if( hypre_CSRMatrixNumRows( A_parcsr->diag ) != hypre_CSRMatrixNumRows( A_parcsr->offd ) )
 {
 MB_SET_ERR_RET( "Row does not match" );
 }
 
 HYPRE_Int size = hypre_CSRMatrixNumRows( A_parcsr->diag );
 HYPRE_Int jj;
 double* Adiag_data = hypre_CSRMatrixData( A_parcsr->diag );
 HYPRE_Int* Adiag_i = hypre_CSRMatrixI( A_parcsr->diag );
 
 for( jj = 0; jj < Adiag_i[size]; ++jj )
 {
 Adiag_data[jj] *= s;
 }
 
 double* Aoffd_data = hypre_CSRMatrixData( A_parcsr->offd );
 HYPRE_Int* Aoffd_i = hypre_CSRMatrixI( A_parcsr->offd );
 
 for( jj = 0; jj < Aoffd_i[size]; ++jj )
 {
 Aoffd_data[jj] *= s;
 }
}
 
HYPRE_Int HypreParMatrix::GetValues( const HYPRE_Int nrows,
 HYPRE_Int* ncols,
 HYPRE_Int* rows,
 HYPRE_Int* cols,
 HYPRE_Complex* values )
{
 return HYPRE_IJMatrixGetValues( A, nrows, ncols, rows, cols, values );
}
 
HYPRE_Int HypreParMatrix::SetValues( const HYPRE_Int nrows,
 HYPRE_Int* ncols,
 const HYPRE_Int* rows,
 const HYPRE_Int* cols,
 const HYPRE_Complex* values )
{
 return HYPRE_IJMatrixSetValues( A, nrows, ncols, rows, cols, values );
}
 
HYPRE_Int HypreParMatrix::AddToValues( const HYPRE_Int nrows,
 HYPRE_Int* ncols,
 const HYPRE_Int* rows,
 const HYPRE_Int* cols,
 const HYPRE_Complex* values )
{
 return HYPRE_IJMatrixAddToValues( A, nrows, ncols, rows, cols, values );
}
 
HYPRE_Int HypreParMatrix::verbosity( const HYPRE_Int level )
{
 return HYPRE_IJMatrixSetPrintLevel( A, level );
}
 
HYPRE_Int HypreParMatrix::FinalizeAssembly()
{
 return HYPRE_IJMatrixAssemble( A );
}
 
static void get_sorted_rows_cols( const std::vector< HYPRE_Int >& rows_cols, std::vector< HYPRE_Int >& hypre_sorted )
{
 hypre_sorted.resize( rows_cols.size() );
 std::copy( rows_cols.begin(), rows_cols.end(), hypre_sorted.begin() );
 std::sort( hypre_sorted.begin(), hypre_sorted.end() );
}
 
void HypreParMatrix::Threshold( double threshold )
{
 int ierr = 0;
 int num_procs;
 hypre_CSRMatrix* csr_A;
 hypre_CSRMatrix* csr_A_wo_z;
 hypre_ParCSRMatrix* parcsr_A_ptr;
 int* row_starts = NULL;
 int* col_starts = NULL;
 int row_start = -1;
 int row_end = -1;
 int col_start = -1;
 int col_end = -1;
 MPI_Comm_size( pcomm->comm(), &num_procs );
 ierr += hypre_ParCSRMatrixGetLocalRange( A_parcsr, &row_start, &row_end, &col_start, &col_end );
 row_starts = hypre_ParCSRMatrixRowStarts( A_parcsr );
 col_starts = hypre_ParCSRMatrixColStarts( A_parcsr );
 parcsr_A_ptr = hypre_ParCSRMatrixCreate( pcomm->comm(), row_starts[num_procs], col_starts[num_procs], row_starts,
 col_starts, 0, 0, 0 );
 csr_A = hypre_MergeDiagAndOffd( A_parcsr );
 csr_A_wo_z = hypre_CSRMatrixDeleteZeros( csr_A, threshold );
 
 /* hypre_CSRMatrixDeleteZeros will return a NULL pointer rather than a usable
 CSR matrix if it finds no non-zeros */
 if( csr_A_wo_z == NULL )
 {
 csr_A_wo_z = csr_A;
 }
 else
 {
 ierr += hypre_CSRMatrixDestroy( csr_A );
 }
 
 ierr += GenerateDiagAndOffd( csr_A_wo_z, parcsr_A_ptr, col_start, col_end );
 ierr += hypre_CSRMatrixDestroy( csr_A_wo_z );
 ierr += hypre_ParCSRMatrixDestroy( A_parcsr );
 hypre_IJMatrixObject( A ) = parcsr_A_ptr;
 /* Get the parcsr matrix object to use */
 HYPRE_IJMatrixGetObject( A, (void**)A_parcsr );
}
 
void HypreParMatrix::EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols,
 const HypreParVector& X,
 HypreParVector& B )
{
 std::vector< HYPRE_Int > rc_sorted;
 get_sorted_rows_cols( rows_cols, rc_sorted );
 internal::hypre_ParCSRMatrixEliminateAXB( A_parcsr, rc_sorted.size(), rc_sorted.data(), X.x_par, B.x_par );
}
 
HypreParMatrix* HypreParMatrix::EliminateRowsCols( const std::vector< HYPRE_Int >& rows_cols )
{
 std::vector< HYPRE_Int > rc_sorted;
 get_sorted_rows_cols( rows_cols, rc_sorted );
 hypre_ParCSRMatrix* Ae;
 internal::hypre_ParCSRMatrixEliminateAAe( A_parcsr, &Ae, rc_sorted.size(), rc_sorted.data() );
 HypreParMatrix* tmpMat = new HypreParMatrix( pcomm, M(), N(), RowPart(), ColPart(), 0, 0 );
 hypre_IJMatrixObject( tmpMat->A ) = Ae;
 /* Get the parcsr matrix object to use */
 HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
 return tmpMat;
}
 
void HypreParMatrix::Print( const char* fname, HYPRE_Int offi, HYPRE_Int offj )
{
 hypre_ParCSRMatrixPrintIJ( A_parcsr, offi, offj, fname );
}
 
void HypreParMatrix::Read( const char* fname )
{
 Destroy();
 Init();
 HYPRE_Int base_i, base_j;
 hypre_ParCSRMatrixReadIJ( pcomm->comm(), fname, &base_i, &base_j, &A_parcsr );
 hypre_ParCSRMatrixSetNumNonzeros( A_parcsr );
 hypre_MatvecCommPkgCreate( A_parcsr );
 height = GetNumRows();
 width = GetNumCols();
}
 
void HypreParMatrix::Destroy()
{
 if( interpreter )
 {
 hypre_TFree( interpreter );
 }
 
 if( A == NULL )
 {
 return;
 }
 
 else
 {
 HYPRE_IJMatrixDestroy( A );
 }
}
 
HypreParMatrix* ParMult( HypreParMatrix* A, HypreParMatrix* B )
{
 hypre_ParCSRMatrix* ab;
 ab = hypre_ParMatmul( A->A_parcsr, B->A_parcsr );
 hypre_MatvecCommPkgCreate( ab );
 HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
 hypre_IJMatrixObject( tmpMat->A ) = ab;
 /* Get the parcsr matrix object to use */
 HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
 return tmpMat;
}
 
HypreParMatrix* RAP( HypreParMatrix* A, HypreParMatrix* P )
{
 HYPRE_Int P_owns_its_col_starts = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( P->A_parcsr ) );
 hypre_ParCSRMatrix* rap;
 hypre_BoomerAMGBuildCoarseOperator( P->A_parcsr, A->A_parcsr, P->A_parcsr, &rap );
 hypre_ParCSRMatrixSetNumNonzeros( rap );
 // hypre_MatvecCommPkgCreate(rap);
 /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
 from P (even if it does not own them)! */
 hypre_ParCSRMatrixSetRowStartsOwner( rap, 0 );
 hypre_ParCSRMatrixSetColStartsOwner( rap, 0 );
 
 if( P_owns_its_col_starts )
 {
 hypre_ParCSRMatrixSetColStartsOwner( P->A_parcsr, 1 );
 }
 
 HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
 hypre_IJMatrixObject( tmpMat->A ) = rap;
 /* Get the parcsr matrix object to use */
 HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
 return tmpMat;
}
 
HypreParMatrix* RAP( HypreParMatrix* Rt, HypreParMatrix* A, HypreParMatrix* P )
{
 HYPRE_Int P_owns_its_col_starts = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( P->A_parcsr ) );
 HYPRE_Int Rt_owns_its_col_starts = hypre_ParCSRMatrixOwnsColStarts( (hypre_ParCSRMatrix*)( Rt->A_parcsr ) );
 hypre_ParCSRMatrix* rap;
 hypre_BoomerAMGBuildCoarseOperator( Rt->A_parcsr, A->A_parcsr, P->A_parcsr, &rap );
 hypre_ParCSRMatrixSetNumNonzeros( rap );
 
 // hypre_MatvecCommPkgCreate(rap);
 if( !P_owns_its_col_starts )
 {
 /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
 from P (even if it does not own them)! */
 hypre_ParCSRMatrixSetColStartsOwner( rap, 0 );
 }
 
 if( !Rt_owns_its_col_starts )
 {
 /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
 from P (even if it does not own them)! */
 hypre_ParCSRMatrixSetRowStartsOwner( rap, 0 );
 }
 
 HypreParMatrix* tmpMat = new HypreParMatrix( A->pcomm, A->M(), A->N(), A->RowPart(), A->ColPart(), 0, 0 );
 hypre_IJMatrixObject( tmpMat->A ) = rap;
 /* Get the parcsr matrix object to use */
 HYPRE_IJMatrixGetObject( tmpMat->A, (void**)tmpMat->A_parcsr );
 return tmpMat;
}
 
void EliminateBC( HypreParMatrix& A,
 HypreParMatrix& Ae,
 const std::vector< int >& ess_dof_list,
 const HypreParVector& X,
 HypreParVector& B )
{
 // B -= Ae*X
 Ae.Mult( -1.0, X, 1.0, B );
 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( (hypre_ParCSRMatrix*)A.A_parcsr );
 double* data = hypre_CSRMatrixData( A_diag );
 HYPRE_Int* I = hypre_CSRMatrixI( A_diag );
#ifdef MOAB_DEBUG
 HYPRE_Int* J = hypre_CSRMatrixJ( A_diag );
 hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( (hypre_ParCSRMatrix*)A.A_parcsr );
 HYPRE_Int* I_offd = hypre_CSRMatrixI( A_offd );
 double* data_offd = hypre_CSRMatrixData( A_offd );
#endif
 std::vector< HYPRE_Complex > bdata( ess_dof_list.size() ), xdata( ess_dof_list.size() );
 B.GetValues( ess_dof_list.size(), ess_dof_list.data(), bdata.data() );
 X.GetValues( ess_dof_list.size(), ess_dof_list.data(), xdata.data() );
 
 for( size_t i = 0; i < ess_dof_list.size(); i++ )
 {
 int r = ess_dof_list[i];
 bdata[r] = data[I[r]] * xdata[r];
#ifdef MOAB_DEBUG
 
 // Check that in the rows specified by the ess_dof_list, the matrix A has
 // only one entry -- the diagonal.
 // if (I[r+1] != I[r]+1 || J[I[r]] != r || I_offd[r] != I_offd[r+1])
 if( J[I[r]] != r )
 {
 MB_SET_ERR_RET( "the diagonal entry must be the first entry in the row!" );
 }
 
 for( int j = I[r] + 1; j < I[r + 1]; j++ )
 {
 if( data[j] != 0.0 )
 {
 MB_SET_ERR_RET( "all off-diagonal entries must be zero!" );
 }
 }
 
 for( int j = I_offd[r]; j < I_offd[r + 1]; j++ )
 {
 if( data_offd[j] != 0.0 )
 {
 MB_SET_ERR_RET( "all off-diagonal entries must be zero!" );
 }
 }
 
#endif
 }
}
 
// Taubin or "lambda-mu" scheme, which alternates between positive and
// negative step sizes to approximate low-pass filter effect.
 
int ParCSRRelax_Taubin( hypre_ParCSRMatrix* A, // matrix to relax with
 hypre_ParVector* f, // right-hand side
 double lambda,
 double mu,
 int N,
 double max_eig,
 hypre_ParVector* u, // initial/updated approximation
 hypre_ParVector* r // another temp vector
)
{
 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
 HYPRE_Int num_rows = hypre_CSRMatrixNumRows( A_diag );
 double* u_data = hypre_VectorData( hypre_ParVectorLocalVector( u ) );
 double* r_data = hypre_VectorData( hypre_ParVectorLocalVector( r ) );
 
 for( int i = 0; i < N; i++ )
 {
 // get residual: r = f - A*u
 hypre_ParVectorCopy( f, r );
 hypre_ParCSRMatrixMatvec( -1.0, A, u, 1.0, r );
 double coef;
 ( 0 == ( i % 2 ) ) ? coef = lambda : coef = mu;
 
 for( HYPRE_Int j = 0; j < num_rows; j++ )
 {
 u_data[j] += coef * r_data[j] / max_eig;
 }
 }
 
 return 0;
}
 
// FIR scheme, which uses Chebyshev polynomials and a window function
// to approximate a low-pass step filter.
 
int ParCSRRelax_FIR( hypre_ParCSRMatrix* A, // matrix to relax with
 hypre_ParVector* f, // right-hand side
 double max_eig,
 int poly_order,
 double* fir_coeffs,
 hypre_ParVector* u, // initial/updated approximation
 hypre_ParVector* x0, // temporaries
 hypre_ParVector* x1,
 hypre_ParVector* x2,
 hypre_ParVector* x3 )
 
{
 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
 HYPRE_Int num_rows = hypre_CSRMatrixNumRows( A_diag );
 double* u_data = hypre_VectorData( hypre_ParVectorLocalVector( u ) );
 double* x0_data = hypre_VectorData( hypre_ParVectorLocalVector( x0 ) );
 double* x1_data = hypre_VectorData( hypre_ParVectorLocalVector( x1 ) );
 double* x2_data = hypre_VectorData( hypre_ParVectorLocalVector( x2 ) );
 double* x3_data = hypre_VectorData( hypre_ParVectorLocalVector( x3 ) );
 hypre_ParVectorCopy( u, x0 );
 // x1 = f -A*x0/max_eig
 hypre_ParVectorCopy( f, x1 );
 hypre_ParCSRMatrixMatvec( -1.0, A, x0, 1.0, x1 );
 
 for( HYPRE_Int i = 0; i < num_rows; i++ )
 {
 x1_data[i] /= -max_eig;
 }
 
 // x1 = x0 -x1
 for( HYPRE_Int i = 0; i < num_rows; i++ )
 {
 x1_data[i] = x0_data[i] - x1_data[i];
 }
 
 // x3 = f0*x0 +f1*x1
 for( HYPRE_Int i = 0; i < num_rows; i++ )
 {
 x3_data[i] = fir_coeffs[0] * x0_data[i] + fir_coeffs[1] * x1_data[i];
 }
 
 for( int n = 2; n <= poly_order; n++ )
 {
 // x2 = f - A*x1/max_eig
 hypre_ParVectorCopy( f, x2 );
 hypre_ParCSRMatrixMatvec( -1.0, A, x1, 1.0, x2 );
 
 for( HYPRE_Int i = 0; i < num_rows; i++ )
 {
 x2_data[i] /= -max_eig;
 }
 
 // x2 = (x1-x0) +(x1-2*x2)
 // x3 = x3 +f[n]*x2
 // x0 = x1
 // x1 = x2
 
 for( HYPRE_Int i = 0; i < num_rows; i++ )
 {
 x2_data[i] = ( x1_data[i] - x0_data[i] ) + ( x1_data[i] - 2 * x2_data[i] );
 x3_data[i] += fir_coeffs[n] * x2_data[i];
 x0_data[i] = x1_data[i];
 x1_data[i] = x2_data[i];
 }
 }
 
 for( HYPRE_Int i = 0; i < num_rows; i++ )
 {
 u_data[i] = x3_data[i];
 }
 
 return 0;
}
 
} // namespace moab
 
#endif