hypre_parcsr.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 "hypre_parcsr.hpp"
 
#ifdef MOAB_HAVE_MPI
 
#include <limits>
 
namespace moab
{
 
namespace internal
{
 
 /*--------------------------------------------------------------------------
 * A*X = B style elimination
 *--------------------------------------------------------------------------*/
 
 /*
 Function: hypre_CSRMatrixEliminateAXB
 
 Eliminate the rows and columns of A corresponding to the
 given sorted (!) list of rows. Put I on the eliminated diagonal,
 subtract columns times X from B.
 */
 void hypre_CSRMatrixEliminateAXB( hypre_CSRMatrix* A,
 HYPRE_Int nrows_to_eliminate,
 HYPRE_Int* rows_to_eliminate,
 hypre_Vector* X,
 hypre_Vector* B )
 {
 HYPRE_Int i, j;
 HYPRE_Int irow, jcol, ibeg, iend, pos;
 HYPRE_Real a;
 
 HYPRE_Int* Ai = hypre_CSRMatrixI( A );
 HYPRE_Int* Aj = hypre_CSRMatrixJ( A );
 HYPRE_Real* Adata = hypre_CSRMatrixData( A );
 HYPRE_Int nrows = hypre_CSRMatrixNumRows( A );
 
 HYPRE_Real* Xdata = hypre_VectorData( X );
 HYPRE_Real* Bdata = hypre_VectorData( B );
 
 /* eliminate the columns */
 for( i = 0; i < nrows; i++ )
 {
 ibeg = Ai[i];
 iend = Ai[i + 1];
 for( j = ibeg; j < iend; j++ )
 {
 jcol = Aj[j];
 pos = hypre_BinarySearch( rows_to_eliminate, jcol, nrows_to_eliminate );
 if( pos != -1 )
 {
 a = Adata[j];
 Adata[j] = 0.0;
 Bdata[i] -= a * Xdata[jcol];
 }
 }
 }
 
 /* remove the rows and set the diagonal equal to 1 */
 for( i = 0; i < nrows_to_eliminate; i++ )
 {
 irow = rows_to_eliminate[i];
 ibeg = Ai[irow];
 iend = Ai[irow + 1];
 for( j = ibeg; j < iend; j++ )
 {
 if( Aj[j] == irow )
 {
 Adata[j] = 1.0;
 }
 else
 {
 Adata[j] = 0.0;
 }
 }
 }
 }
 
 /*
 Function: hypre_CSRMatrixEliminateOffdColsAXB
 
 Eliminate the given sorted (!) list of columns of A, subtract them from B.
 */
 void hypre_CSRMatrixEliminateOffdColsAXB( hypre_CSRMatrix* A,
 HYPRE_Int ncols_to_eliminate,
 HYPRE_Int* eliminate_cols,
 HYPRE_Real* eliminate_coefs,
 hypre_Vector* B )
 {
 HYPRE_Int i, j;
 HYPRE_Int ibeg, iend, pos;
 HYPRE_Real a;
 
 HYPRE_Int* Ai = hypre_CSRMatrixI( A );
 HYPRE_Int* Aj = hypre_CSRMatrixJ( A );
 HYPRE_Real* Adata = hypre_CSRMatrixData( A );
 HYPRE_Int nrows = hypre_CSRMatrixNumRows( A );
 
 HYPRE_Real* Bdata = hypre_VectorData( B );
 
 for( i = 0; i < nrows; i++ )
 {
 ibeg = Ai[i];
 iend = Ai[i + 1];
 for( j = ibeg; j < iend; j++ )
 {
 pos = hypre_BinarySearch( eliminate_cols, Aj[j], ncols_to_eliminate );
 if( pos != -1 )
 {
 a = Adata[j];
 Adata[j] = 0.0;
 Bdata[i] -= a * eliminate_coefs[pos];
 }
 }
 }
 }
 
 /*
 Function: hypre_CSRMatrixEliminateOffdRowsAXB
 
 Eliminate (zero) the given list of rows of A.
 */
 void hypre_CSRMatrixEliminateOffdRowsAXB( hypre_CSRMatrix* A,
 HYPRE_Int nrows_to_eliminate,
 HYPRE_Int* rows_to_eliminate )
 {
 HYPRE_Int* Ai = hypre_CSRMatrixI( A );
 HYPRE_Real* Adata = hypre_CSRMatrixData( A );
 
 HYPRE_Int i, j;
 HYPRE_Int irow, ibeg, iend;
 
 for( i = 0; i < nrows_to_eliminate; i++ )
 {
 irow = rows_to_eliminate[i];
 ibeg = Ai[irow];
 iend = Ai[irow + 1];
 for( j = ibeg; j < iend; j++ )
 {
 Adata[j] = 0.0;
 }
 }
 }
 
 /*
 Function: hypre_ParCSRMatrixEliminateAXB
 
 This function eliminates the global rows and columns of a matrix
 A corresponding to given lists of sorted (!) local row numbers,
 so that the solution to the system A*X = B is X_b for the given rows.
 
 The elimination is done as follows:
 
 (input) (output)
 
 / A_ii | A_ib \ / A_ii | 0 \
 A = | -----+----- | ---> | -----+----- |
 \ A_bi | A_bb / \ 0 | I /
 
 / X_i \ / X_i \
 X = | --- | ---> | --- | (no change)
 \ X_b / \ X_b /
 
 / B_i \ / B_i - A_ib * X_b \
 B = | --- | ---> | ---------------- |
 \ B_b / \ X_b /
 
 */
 void hypre_ParCSRMatrixEliminateAXB( hypre_ParCSRMatrix* A,
 HYPRE_Int num_rowscols_to_elim,
 HYPRE_Int* rowscols_to_elim,
 hypre_ParVector* X,
 hypre_ParVector* B )
 {
 hypre_CSRMatrix* diag = hypre_ParCSRMatrixDiag( A );
 hypre_CSRMatrix* offd = hypre_ParCSRMatrixOffd( A );
 HYPRE_Int diag_nrows = hypre_CSRMatrixNumRows( diag );
 HYPRE_Int offd_ncols = hypre_CSRMatrixNumCols( offd );
 
 hypre_Vector* Xlocal = hypre_ParVectorLocalVector( X );
 hypre_Vector* Blocal = hypre_ParVectorLocalVector( B );
 
 HYPRE_Real* Bdata = hypre_VectorData( Blocal );
 HYPRE_Real* Xdata = hypre_VectorData( Xlocal );
 
 HYPRE_Int num_offd_cols_to_elim;
 HYPRE_Int* offd_cols_to_elim;
 HYPRE_Real* eliminate_coefs;
 
 /* figure out which offd cols should be eliminated and with what coef */
 hypre_ParCSRCommHandle* comm_handle;
 hypre_ParCSRCommPkg* comm_pkg;
 HYPRE_Int num_sends;
 HYPRE_Int index, start;
 HYPRE_Int i, j, k, irow;
 
 HYPRE_Real* eliminate_row = hypre_CTAlloc( HYPRE_Real, diag_nrows );
 HYPRE_Real* eliminate_col = hypre_CTAlloc( HYPRE_Real, offd_ncols );
 HYPRE_Real *buf_data, coef;
 
 /* make sure A has a communication package */
 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
 if( !comm_pkg )
 {
 hypre_MatvecCommPkgCreate( A );
 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
 }
 
 /* HACK: rows that shouldn't be eliminated are marked with quiet NaN;
 those that should are set to the boundary value from X; this is to
 avoid sending complex type (int+double) or communicating twice. */
 for( i = 0; i < diag_nrows; i++ )
 {
 eliminate_row[i] = std::numeric_limits< HYPRE_Real >::quiet_NaN();
 }
 for( i = 0; i < num_rowscols_to_elim; i++ )
 {
 irow = rowscols_to_elim[i];
 eliminate_row[irow] = Xdata[irow];
 }
 
 /* use a Matvec communication pattern to find (in eliminate_col)
 which of the local offd columns are to be eliminated */
 num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
 buf_data = hypre_CTAlloc( HYPRE_Real, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
 index = 0;
 for( i = 0; i < num_sends; i++ )
 {
 start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
 for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
 {
 k = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
 buf_data[index++] = eliminate_row[k];
 }
 }
 comm_handle = hypre_ParCSRCommHandleCreate( 1, comm_pkg, buf_data, eliminate_col );
 
 /* do sequential part of the elimination while stuff is getting sent */
 hypre_CSRMatrixEliminateAXB( diag, num_rowscols_to_elim, rowscols_to_elim, Xlocal, Blocal );
 
 /* finish the communication */
 hypre_ParCSRCommHandleDestroy( comm_handle );
 
 /* received eliminate_col[], count offd columns to eliminate */
 num_offd_cols_to_elim = 0;
 for( i = 0; i < offd_ncols; i++ )
 {
 coef = eliminate_col[i];
 if( coef == coef ) // test for NaN
 {
 num_offd_cols_to_elim++;
 }
 }
 
 offd_cols_to_elim = hypre_CTAlloc( HYPRE_Int, num_offd_cols_to_elim );
 eliminate_coefs = hypre_CTAlloc( HYPRE_Real, num_offd_cols_to_elim );
 
 /* get a list of offd column indices and coefs */
 num_offd_cols_to_elim = 0;
 for( i = 0; i < offd_ncols; i++ )
 {
 coef = eliminate_col[i];
 if( coef == coef ) // test for NaN
 {
 offd_cols_to_elim[num_offd_cols_to_elim] = i;
 eliminate_coefs[num_offd_cols_to_elim] = coef;
 num_offd_cols_to_elim++;
 }
 }
 
 hypre_TFree( buf_data );
 hypre_TFree( eliminate_row );
 hypre_TFree( eliminate_col );
 
 /* eliminate the off-diagonal part */
 hypre_CSRMatrixEliminateOffdColsAXB( offd, num_offd_cols_to_elim, offd_cols_to_elim, eliminate_coefs, Blocal );
 
 hypre_CSRMatrixEliminateOffdRowsAXB( offd, num_rowscols_to_elim, rowscols_to_elim );
 
 /* set boundary values in the rhs */
 for( int m = 0; m < num_rowscols_to_elim; m++ )
 {
 irow = rowscols_to_elim[m];
 Bdata[irow] = Xdata[irow];
 }
 
 hypre_TFree( offd_cols_to_elim );
 hypre_TFree( eliminate_coefs );
 }
 
 /*--------------------------------------------------------------------------
 * (A + Ae) style elimination
 *--------------------------------------------------------------------------*/
 
 /*
 Function: hypre_CSRMatrixElimCreate
 
 Prepare the Ae matrix: count nnz, initialize I, allocate J and data.
 */
 void hypre_CSRMatrixElimCreate( hypre_CSRMatrix* A,
 hypre_CSRMatrix* Ae,
 HYPRE_Int nrows,
 HYPRE_Int* rows,
 HYPRE_Int ncols,
 HYPRE_Int* cols,
 HYPRE_Int* col_mark )
 {
 HYPRE_Int i, j, col;
 HYPRE_Int A_beg, A_end;
 
 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
 HYPRE_Int A_rows = hypre_CSRMatrixNumRows( A );
 
 hypre_CSRMatrixI( Ae ) = hypre_TAlloc( HYPRE_Int, A_rows + 1 );
 
 HYPRE_Int* Ae_i = hypre_CSRMatrixI( Ae );
 HYPRE_Int nnz = 0;
 
 for( i = 0; i < A_rows; i++ )
 {
 Ae_i[i] = nnz;
 
 A_beg = A_i[i];
 A_end = A_i[i + 1];
 
 if( hypre_BinarySearch( rows, i, nrows ) >= 0 )
 {
 /* full row */
 nnz += A_end - A_beg;
 
 if( col_mark )
 {
 for( j = A_beg; j < A_end; j++ )
 {
 col_mark[A_j[j]] = 1;
 }
 }
 }
 else
 {
 /* count columns */
 for( j = A_beg; j < A_end; j++ )
 {
 col = A_j[j];
 if( hypre_BinarySearch( cols, col, ncols ) >= 0 )
 {
 nnz++;
 if( col_mark )
 {
 col_mark[col] = 1;
 }
 }
 }
 }
 }
 Ae_i[A_rows] = nnz;
 
 hypre_CSRMatrixJ( Ae ) = hypre_TAlloc( HYPRE_Int, nnz );
 hypre_CSRMatrixData( Ae ) = hypre_TAlloc( HYPRE_Real, nnz );
 hypre_CSRMatrixNumNonzeros( Ae ) = nnz;
 }
 
 /*
 Function: hypre_CSRMatrixEliminateRowsCols
 
 Eliminate rows and columns of A, store eliminated values in Ae.
 If 'diag' is nonzero, the eliminated diagonal of A is set to identity.
 If 'col_remap' is not NULL it specifies renumbering of columns of Ae.
 */
 void hypre_CSRMatrixEliminateRowsCols( hypre_CSRMatrix* A,
 hypre_CSRMatrix* Ae,
 HYPRE_Int nrows,
 HYPRE_Int* rows,
 HYPRE_Int ncols,
 HYPRE_Int* cols,
 int diag,
 HYPRE_Int* col_remap )
 {
 HYPRE_Int i, j, k, col;
 HYPRE_Int A_beg, Ae_beg, A_end;
 HYPRE_Real a;
 
 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
 HYPRE_Real* A_data = hypre_CSRMatrixData( A );
 HYPRE_Int A_rows = hypre_CSRMatrixNumRows( A );
 
 HYPRE_Int* Ae_i = hypre_CSRMatrixI( Ae );
 HYPRE_Int* Ae_j = hypre_CSRMatrixJ( Ae );
 HYPRE_Real* Ae_data = hypre_CSRMatrixData( Ae );
 
 for( i = 0; i < A_rows; i++ )
 {
 A_beg = A_i[i];
 A_end = A_i[i + 1];
 Ae_beg = Ae_i[i];
 
 if( hypre_BinarySearch( rows, i, nrows ) >= 0 )
 {
 /* eliminate row */
 for( j = A_beg, k = Ae_beg; j < A_end; j++, k++ )
 {
 col = A_j[j];
 Ae_j[k] = col_remap ? col_remap[col] : col;
 a = ( diag && col == i ) ? 1.0 : 0.0;
 Ae_data[k] = A_data[j] - a;
 A_data[j] = a;
 }
 }
 else
 {
 /* eliminate columns */
 for( j = A_beg, k = Ae_beg; j < A_end; j++ )
 {
 col = A_j[j];
 if( hypre_BinarySearch( cols, col, ncols ) >= 0 )
 {
 Ae_j[k] = col_remap ? col_remap[col] : col;
 Ae_data[k] = A_data[j];
 A_data[j] = 0.0;
 k++;
 }
 }
 }
 }
 }
 
 /*
 Function: hypre_ParCSRMatrixEliminateAAe
 
 (input) (output)
 
 / A_ii | A_ib \ / A_ii | 0 \
 A = | -----+----- | ---> | -----+----- |
 \ A_bi | A_bb / \ 0 | I /
 
 
 / 0 | A_ib
 \ Ae = | -----+--------- | \ A_bi | A_bb - I /
 
 */
 void hypre_ParCSRMatrixEliminateAAe( hypre_ParCSRMatrix* A,
 hypre_ParCSRMatrix** Ae,
 HYPRE_Int num_rowscols_to_elim,
 HYPRE_Int* rowscols_to_elim )
 {
 HYPRE_Int i, j, k;
 
 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
 hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
 HYPRE_Int A_diag_nrows = hypre_CSRMatrixNumRows( A_diag );
 HYPRE_Int A_offd_ncols = hypre_CSRMatrixNumCols( A_offd );
 
 *Ae = hypre_ParCSRMatrixCreate( hypre_ParCSRMatrixComm( A ), hypre_ParCSRMatrixGlobalNumRows( A ),
 hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
 hypre_ParCSRMatrixColStarts( A ), 0, 0, 0 );
 
 hypre_ParCSRMatrixSetRowStartsOwner( *Ae, 0 );
 hypre_ParCSRMatrixSetColStartsOwner( *Ae, 0 );
 
 hypre_CSRMatrix* Ae_diag = hypre_ParCSRMatrixDiag( *Ae );
 hypre_CSRMatrix* Ae_offd = hypre_ParCSRMatrixOffd( *Ae );
 HYPRE_Int Ae_offd_ncols;
 
 HYPRE_Int num_offd_cols_to_elim;
 HYPRE_Int* offd_cols_to_elim;
 
 HYPRE_Int* A_col_map_offd = hypre_ParCSRMatrixColMapOffd( A );
 HYPRE_Int* Ae_col_map_offd;
 
 HYPRE_Int* col_mark;
 HYPRE_Int* col_remap;
 
 /* figure out which offd cols should be eliminated */
 {
 hypre_ParCSRCommHandle* comm_handle;
 hypre_ParCSRCommPkg* comm_pkg;
 HYPRE_Int num_sends, *int_buf_data;
 HYPRE_Int index, start;
 
 HYPRE_Int* eliminate_row = hypre_CTAlloc( HYPRE_Int, A_diag_nrows );
 HYPRE_Int* eliminate_col = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
 
 /* make sure A has a communication package */
 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
 if( !comm_pkg )
 {
 hypre_MatvecCommPkgCreate( A );
 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
 }
 
 /* which of the local rows are to be eliminated */
 for( i = 0; i < A_diag_nrows; i++ )
 {
 eliminate_row[i] = 0;
 }
 for( i = 0; i < num_rowscols_to_elim; i++ )
 {
 eliminate_row[rowscols_to_elim[i]] = 1;
 }
 
 /* use a Matvec communication pattern to find (in eliminate_col)
 which of the local offd columns are to be eliminated */
 num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
 int_buf_data = hypre_CTAlloc( HYPRE_Int, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
 index = 0;
 for( i = 0; i < num_sends; i++ )
 {
 start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
 for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
 {
 k = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
 int_buf_data[index++] = eliminate_row[k];
 }
 }
 comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data, eliminate_col );
 
 /* eliminate diagonal part, overlapping it with communication */
 hypre_CSRMatrixElimCreate( A_diag, Ae_diag, num_rowscols_to_elim, rowscols_to_elim, num_rowscols_to_elim,
 rowscols_to_elim, NULL );
 
 hypre_CSRMatrixEliminateRowsCols( A_diag, Ae_diag, num_rowscols_to_elim, rowscols_to_elim,
 num_rowscols_to_elim, rowscols_to_elim, 1, NULL );
 hypre_CSRMatrixReorder( Ae_diag );
 
 /* finish the communication */
 hypre_ParCSRCommHandleDestroy( comm_handle );
 
 /* received eliminate_col[], count offd columns to eliminate */
 num_offd_cols_to_elim = 0;
 for( i = 0; i < A_offd_ncols; i++ )
 {
 if( eliminate_col[i] )
 {
 num_offd_cols_to_elim++;
 }
 }
 
 offd_cols_to_elim = hypre_CTAlloc( HYPRE_Int, num_offd_cols_to_elim );
 
 /* get a list of offd column indices and coefs */
 num_offd_cols_to_elim = 0;
 for( i = 0; i < A_offd_ncols; i++ )
 {
 if( eliminate_col[i] )
 {
 offd_cols_to_elim[num_offd_cols_to_elim++] = i;
 }
 }
 
 hypre_TFree( int_buf_data );
 hypre_TFree( eliminate_row );
 hypre_TFree( eliminate_col );
 }
 
 /* eliminate the off-diagonal part */
 col_mark = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
 col_remap = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
 
 hypre_CSRMatrixElimCreate( A_offd, Ae_offd, num_rowscols_to_elim, rowscols_to_elim, num_offd_cols_to_elim,
 offd_cols_to_elim, col_mark );
 
 for( i = k = 0; i < A_offd_ncols; i++ )
 {
 if( col_mark[i] )
 {
 col_remap[i] = k++;
 }
 }
 
 hypre_CSRMatrixEliminateRowsCols( A_offd, Ae_offd, num_rowscols_to_elim, rowscols_to_elim,
 num_offd_cols_to_elim, offd_cols_to_elim, 0, col_remap );
 
 /* create col_map_offd for Ae */
 Ae_offd_ncols = 0;
 for( i = 0; i < A_offd_ncols; i++ )
 {
 if( col_mark[i] )
 {
 Ae_offd_ncols++;
 }
 }
 
 Ae_col_map_offd = hypre_CTAlloc( HYPRE_Int, Ae_offd_ncols );
 
 Ae_offd_ncols = 0;
 for( i = 0; i < A_offd_ncols; i++ )
 {
 if( col_mark[i] )
 {
 Ae_col_map_offd[Ae_offd_ncols++] = A_col_map_offd[i];
 }
 }
 
 hypre_ParCSRMatrixColMapOffd( *Ae ) = Ae_col_map_offd;
 hypre_CSRMatrixNumCols( Ae_offd ) = Ae_offd_ncols;
 
 hypre_TFree( col_remap );
 hypre_TFree( col_mark );
 hypre_TFree( offd_cols_to_elim );
 
 hypre_ParCSRMatrixSetNumNonzeros( *Ae );
 hypre_MatvecCommPkgCreate( *Ae );
 }
 
 /*--------------------------------------------------------------------------
 * Split
 *--------------------------------------------------------------------------*/
 
 void hypre_CSRMatrixSplit( hypre_CSRMatrix* A,
 HYPRE_Int nr,
 HYPRE_Int nc,
 HYPRE_Int* row_block_num,
 HYPRE_Int* col_block_num,
 hypre_CSRMatrix** blocks )
 {
 HYPRE_Int i, j, k, bi, bj;
 
 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
 HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
 
 HYPRE_Int A_rows = hypre_CSRMatrixNumRows( A );
 HYPRE_Int A_cols = hypre_CSRMatrixNumCols( A );
 
 HYPRE_Int* num_rows = hypre_CTAlloc( HYPRE_Int, nr );
 HYPRE_Int* num_cols = hypre_CTAlloc( HYPRE_Int, nc );
 
 HYPRE_Int* block_row = hypre_TAlloc( HYPRE_Int, A_rows );
 HYPRE_Int* block_col = hypre_TAlloc( HYPRE_Int, A_cols );
 
 for( i = 0; i < A_rows; i++ )
 {
 block_row[i] = num_rows[row_block_num[i]]++;
 }
 for( j = 0; j < A_cols; j++ )
 {
 block_col[j] = num_cols[col_block_num[j]]++;
 }
 
 /* allocate the blocks */
 for( i = 0; i < nr; i++ )
 {
 for( j = 0; j < nc; j++ )
 {
 hypre_CSRMatrix* B = hypre_CSRMatrixCreate( num_rows[i], num_cols[j], 0 );
 hypre_CSRMatrixI( B ) = hypre_CTAlloc( HYPRE_Int, num_rows[i] + 1 );
 blocks[i * nc + j] = B;
 }
 }
 
 /* count block row nnz */
 for( i = 0; i < A_rows; i++ )
 {
 bi = row_block_num[i];
 for( j = A_i[i]; j < A_i[i + 1]; j++ )
 {
 bj = col_block_num[A_j[j]];
 hypre_CSRMatrix* B = blocks[bi * nc + bj];
 hypre_CSRMatrixI( B )[block_row[i] + 1]++;
 }
 }
 
 /* count block nnz */
 for( k = 0; k < nr * nc; k++ )
 {
 hypre_CSRMatrix* B = blocks[k];
 HYPRE_Int* B_i = hypre_CSRMatrixI( B );
 
 HYPRE_Int nnz = 0, rs;
 for( int m = 1; m <= hypre_CSRMatrixNumRows( B ); m++ )
 {
 rs = B_i[m], B_i[m] = nnz, nnz += rs;
 }
 
 hypre_CSRMatrixJ( B ) = hypre_TAlloc( HYPRE_Int, nnz );
 hypre_CSRMatrixData( B ) = hypre_TAlloc( HYPRE_Complex, nnz );
 hypre_CSRMatrixNumNonzeros( B ) = nnz;
 }
 
 /* populate blocks */
 for( i = 0; i < A_rows; i++ )
 {
 bi = row_block_num[i];
 for( j = A_i[i]; j < A_i[i + 1]; j++ )
 {
 k = A_j[j];
 bj = col_block_num[k];
 hypre_CSRMatrix* B = blocks[bi * nc + bj];
 HYPRE_Int* bii = hypre_CSRMatrixI( B ) + block_row[i] + 1;
 hypre_CSRMatrixJ( B )[*bii] = block_col[k];
 hypre_CSRMatrixData( B )[*bii] = A_data[j];
 ( *bii )++;
 }
 }
 
 hypre_TFree( block_col );
 hypre_TFree( block_row );
 
 hypre_TFree( num_cols );
 hypre_TFree( num_rows );
 }
 
 void hypre_ParCSRMatrixSplit( hypre_ParCSRMatrix* A,
 HYPRE_Int nr,
 HYPRE_Int nc,
 hypre_ParCSRMatrix** blocks,
 int interleaved_rows,
 int interleaved_cols )
 {
 HYPRE_Int i, j, k;
 
 MPI_Comm comm = hypre_ParCSRMatrixComm( A );
 
 hypre_CSRMatrix* Adiag = hypre_ParCSRMatrixDiag( A );
 hypre_CSRMatrix* Aoffd = hypre_ParCSRMatrixOffd( A );
 
 HYPRE_Int global_rows = hypre_ParCSRMatrixGlobalNumRows( A );
 HYPRE_Int global_cols = hypre_ParCSRMatrixGlobalNumCols( A );
 
 HYPRE_Int local_rows = hypre_CSRMatrixNumRows( Adiag );
 HYPRE_Int local_cols = hypre_CSRMatrixNumCols( Adiag );
 HYPRE_Int offd_cols = hypre_CSRMatrixNumCols( Aoffd );
 
 hypre_assert( local_rows % nr == 0 && local_cols % nc == 0 );
 hypre_assert( global_rows % nr == 0 && global_cols % nc == 0 );
 
 HYPRE_Int block_rows = local_rows / nr;
 HYPRE_Int block_cols = local_cols / nc;
 HYPRE_Int num_blocks = nr * nc;
 
 /* mark local rows and columns with block number */
 HYPRE_Int* row_block_num = hypre_TAlloc( HYPRE_Int, local_rows );
 HYPRE_Int* col_block_num = hypre_TAlloc( HYPRE_Int, local_cols );
 
 for( i = 0; i < local_rows; i++ )
 {
 row_block_num[i] = interleaved_rows ? ( i % nr ) : ( i / block_rows );
 }
 for( i = 0; i < local_cols; i++ )
 {
 col_block_num[i] = interleaved_cols ? ( i % nc ) : ( i / block_cols );
 }
 
 /* determine the block numbers for offd columns */
 HYPRE_Int* offd_col_block_num = hypre_TAlloc( HYPRE_Int, offd_cols );
 hypre_ParCSRCommHandle* comm_handle;
 HYPRE_Int* int_buf_data;
 {
 /* make sure A has a communication package */
 hypre_ParCSRCommPkg* comm_pkg = hypre_ParCSRMatrixCommPkg( A );
 if( !comm_pkg )
 {
 hypre_MatvecCommPkgCreate( A );
 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
 }
 
 /* calculate the final global column numbers for each block */
 HYPRE_Int* count = hypre_CTAlloc( HYPRE_Int, nc );
 HYPRE_Int* block_global_col = hypre_TAlloc( HYPRE_Int, local_cols );
 HYPRE_Int first_col = hypre_ParCSRMatrixFirstColDiag( A ) / nc;
 for( i = 0; i < local_cols; i++ )
 {
 block_global_col[i] = first_col + count[col_block_num[i]]++;
 }
 hypre_TFree( count );
 
 /* use a Matvec communication pattern to determine offd_col_block_num */
 HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
 int_buf_data = hypre_CTAlloc( HYPRE_Int, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
 HYPRE_Int start, index = 0;
 for( i = 0; i < num_sends; i++ )
 {
 start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
 for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
 {
 k = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
 int_buf_data[index++] = col_block_num[k] + nc * block_global_col[k];
 }
 }
 hypre_TFree( block_global_col );
 
 comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data, offd_col_block_num );
 }
 
 /* create the block matrices */
 HYPRE_Int num_procs = 1;
 if( !hypre_ParCSRMatrixAssumedPartition( A ) )
 {
 hypre_MPI_Comm_size( comm, &num_procs );
 }
 
 HYPRE_Int* row_starts = hypre_TAlloc( HYPRE_Int, num_procs + 1 );
 HYPRE_Int* col_starts = hypre_TAlloc( HYPRE_Int, num_procs + 1 );
 for( i = 0; i <= num_procs; i++ )
 {
 row_starts[i] = hypre_ParCSRMatrixRowStarts( A )[i] / nr;
 col_starts[i] = hypre_ParCSRMatrixColStarts( A )[i] / nc;
 }
 
 for( i = 0; i < num_blocks; i++ )
 {
 blocks[i] =
 hypre_ParCSRMatrixCreate( comm, global_rows / nr, global_cols / nc, row_starts, col_starts, 0, 0, 0 );
 }
 
 /* split diag part */
 hypre_CSRMatrix** csr_blocks = hypre_TAlloc( hypre_CSRMatrix*, nr * nc );
 hypre_CSRMatrixSplit( Adiag, nr, nc, row_block_num, col_block_num, csr_blocks );
 
 for( i = 0; i < num_blocks; i++ )
 {
 hypre_TFree( hypre_ParCSRMatrixDiag( blocks[i] ) );
 hypre_ParCSRMatrixDiag( blocks[i] ) = csr_blocks[i];
 }
 
 /* finish communication, receive offd_col_block_num */
 hypre_ParCSRCommHandleDestroy( comm_handle );
 hypre_TFree( int_buf_data );
 
 /* decode global offd column numbers */
 HYPRE_Int* offd_global_col = hypre_TAlloc( HYPRE_Int, offd_cols );
 for( i = 0; i < offd_cols; i++ )
 {
 offd_global_col[i] = offd_col_block_num[i] / nc;
 offd_col_block_num[i] %= nc;
 }
 
 /* split offd part */
 hypre_CSRMatrixSplit( Aoffd, nr, nc, row_block_num, offd_col_block_num, csr_blocks );
 
 for( i = 0; i < num_blocks; i++ )
 {
 hypre_TFree( hypre_ParCSRMatrixOffd( blocks[i] ) );
 hypre_ParCSRMatrixOffd( blocks[i] ) = csr_blocks[i];
 }
 
 hypre_TFree( csr_blocks );
 hypre_TFree( col_block_num );
 hypre_TFree( row_block_num );
 
 /* update block col-maps */
 for( int bi = 0; bi < nr; bi++ )
 {
 for( int bj = 0; bj < nc; bj++ )
 {
 hypre_ParCSRMatrix* block = blocks[bi * nc + bj];
 hypre_CSRMatrix* block_offd = hypre_ParCSRMatrixOffd( block );
 HYPRE_Int block_offd_cols = hypre_CSRMatrixNumCols( block_offd );
 
 HYPRE_Int* block_col_map = hypre_TAlloc( HYPRE_Int, block_offd_cols );
 for( i = j = 0; i < offd_cols; i++ )
 {
 HYPRE_Int bn = offd_col_block_num[i];
 if( bn == bj )
 {
 block_col_map[j++] = offd_global_col[i];
 }
 }
 hypre_assert( j == block_offd_cols );
 
 hypre_ParCSRMatrixColMapOffd( block ) = block_col_map;
 }
 }
 
 hypre_TFree( offd_global_col );
 hypre_TFree( offd_col_block_num );
 
 /* finish the new matrices, make them own all the stuff */
 for( i = 0; i < num_blocks; i++ )
 {
 hypre_ParCSRMatrixSetNumNonzeros( blocks[i] );
 hypre_MatvecCommPkgCreate( blocks[i] );
 
 hypre_ParCSRMatrixOwnsData( blocks[i] ) = 1;
 
 /* only the first block will own the row/col_starts */
 hypre_ParCSRMatrixOwnsRowStarts( blocks[i] ) = !i;
 hypre_ParCSRMatrixOwnsColStarts( blocks[i] ) = !i;
 }
 }
 
 /* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
 void hypre_CSRMatrixBooleanMatvec( hypre_CSRMatrix* A,
 HYPRE_Bool alpha,
 HYPRE_Bool* x,
 HYPRE_Bool beta,
 HYPRE_Bool* y )
 {
 /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
 HYPRE_Int num_rows = hypre_CSRMatrixNumRows( A );
 
 HYPRE_Int* A_rownnz = hypre_CSRMatrixRownnz( A );
 HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz( A );
 
 HYPRE_Bool* x_data = x;
 HYPRE_Bool* y_data = y;
 
 HYPRE_Bool temp, tempx;
 
 HYPRE_Int i, jj;
 
 HYPRE_Int m;
 
 HYPRE_Real xpar = 0.7;
 
 /*-----------------------------------------------------------------------
 * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
 *-----------------------------------------------------------------------*/
 
 if( alpha == 0 )
 {
#ifdef HYPRE_USING_OPENMP
#pragma omp parallel for private( i ) HYPRE_SMP_SCHEDULE
#endif
 for( i = 0; i < num_rows; i++ )
 {
 y_data[i] = y_data[i] && beta;
 }
 return;
 }
 
 /*-----------------------------------------------------------------------
 * y = (beta/alpha)*y
 *-----------------------------------------------------------------------*/
 
 if( beta == 0 )
 {
#ifdef HYPRE_USING_OPENMP
#pragma omp parallel for private( i ) HYPRE_SMP_SCHEDULE
#endif
 for( i = 0; i < num_rows; i++ )
 {
 y_data[i] = 0;
 }
 }
 else
 {
 /* beta is true -> no change to y_data */
 }
 
 /*-----------------------------------------------------------------
 * y += A*x
 *-----------------------------------------------------------------*/
 
 /* use rownnz pointer to do the A*x multiplication when num_rownnz is smaller than num_rows
 */
 
 if( num_rownnz < xpar * ( num_rows ) )
 {
#ifdef HYPRE_USING_OPENMP
#pragma omp parallel for private( i, jj, m, tempx ) HYPRE_SMP_SCHEDULE
#endif
 for( i = 0; i < num_rownnz; i++ )
 {
 m = A_rownnz[i];
 
 tempx = 0;
 for( jj = A_i[m]; jj < A_i[m + 1]; jj++ )
 {
 /* tempx = tempx || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
 tempx = tempx || x_data[A_j[jj]];
 }
 y_data[m] = y_data[m] || tempx;
 }
 }
 else
 {
#ifdef HYPRE_USING_OPENMP
#pragma omp parallel for private( i, jj, temp ) HYPRE_SMP_SCHEDULE
#endif
 for( i = 0; i < num_rows; i++ )
 {
 temp = 0;
 for( jj = A_i[i]; jj < A_i[i + 1]; jj++ )
 {
 /* temp = temp || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
 temp = temp || x_data[A_j[jj]];
 }
 y_data[i] = y_data[i] || temp;
 }
 }
 
 /*-----------------------------------------------------------------
 * y = alpha*y
 *-----------------------------------------------------------------*/
 /* alpha is true */
 }
 
 /* Based on hypre_ParCSRCommHandleCreate in hypre's par_csr_communication.c. The
 input variable job controls the communication type: 1=Matvec, 2=MatvecT. */
 hypre_ParCSRCommHandle* hypre_ParCSRCommHandleCreate_bool( HYPRE_Int job,
 hypre_ParCSRCommPkg* comm_pkg,
 HYPRE_Bool* send_data,
 HYPRE_Bool* recv_data )
 {
 HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
 HYPRE_Int num_recvs = hypre_ParCSRCommPkgNumRecvs( comm_pkg );
 MPI_Comm comm = hypre_ParCSRCommPkgComm( comm_pkg );
 
 hypre_ParCSRCommHandle* comm_handle;
 HYPRE_Int num_requests;
 hypre_MPI_Request* requests;
 
 HYPRE_Int i, j;
 HYPRE_Int my_id, num_procs;
 HYPRE_Int ip, vec_start, vec_len;
 
 num_requests = num_sends + num_recvs;
 requests = hypre_CTAlloc( hypre_MPI_Request, num_requests );
 
 hypre_MPI_Comm_size( comm, &num_procs );
 hypre_MPI_Comm_rank( comm, &my_id );
 
 j = 0;
 switch( job )
 {
 case 1: {
 HYPRE_Bool* d_send_data = (HYPRE_Bool*)send_data;
 HYPRE_Bool* d_recv_data = (HYPRE_Bool*)recv_data;
 for( i = 0; i < num_recvs; i++ )
 {
 ip = hypre_ParCSRCommPkgRecvProc( comm_pkg, i );
 vec_start = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i );
 vec_len = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i + 1 ) - vec_start;
 hypre_MPI_Irecv( &d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
 }
 for( i = 0; i < num_sends; i++ )
 {
 vec_start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
 vec_len = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ) - vec_start;
 ip = hypre_ParCSRCommPkgSendProc( comm_pkg, i );
 hypre_MPI_Isend( &d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
 }
 break;
 }
 case 2: {
 HYPRE_Bool* d_send_data = (HYPRE_Bool*)send_data;
 HYPRE_Bool* d_recv_data = (HYPRE_Bool*)recv_data;
 for( i = 0; i < num_sends; i++ )
 {
 vec_start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
 vec_len = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ) - vec_start;
 ip = hypre_ParCSRCommPkgSendProc( comm_pkg, i );
 hypre_MPI_Irecv( &d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
 }
 for( i = 0; i < num_recvs; i++ )
 {
 ip = hypre_ParCSRCommPkgRecvProc( comm_pkg, i );
 vec_start = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i );
 vec_len = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i + 1 ) - vec_start;
 hypre_MPI_Isend( &d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
 }
 break;
 }
 }
 /*--------------------------------------------------------------------
 * set up comm_handle and return
 *--------------------------------------------------------------------*/
 
 comm_handle = hypre_CTAlloc( hypre_ParCSRCommHandle, 1 );
 
 hypre_ParCSRCommHandleCommPkg( comm_handle ) = comm_pkg;
 hypre_ParCSRCommHandleSendData( comm_handle ) = send_data;
 hypre_ParCSRCommHandleRecvData( comm_handle ) = recv_data;
 hypre_ParCSRCommHandleNumRequests( comm_handle ) = num_requests;
 hypre_ParCSRCommHandleRequests( comm_handle ) = requests;
 
 return comm_handle;
 }
 
 /* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
 void hypre_ParCSRMatrixBooleanMatvec( hypre_ParCSRMatrix* A,
 HYPRE_Bool alpha,
 HYPRE_Bool* x,
 HYPRE_Bool beta,
 HYPRE_Bool* y )
 {
 hypre_ParCSRCommHandle* comm_handle;
 hypre_ParCSRCommPkg* comm_pkg = hypre_ParCSRMatrixCommPkg( A );
 hypre_CSRMatrix* diag = hypre_ParCSRMatrixDiag( A );
 hypre_CSRMatrix* offd = hypre_ParCSRMatrixOffd( A );
 
 HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols( offd );
 HYPRE_Int num_sends, i, j, index;
 
 HYPRE_Bool *x_tmp, *x_buf;
 
 x_tmp = hypre_CTAlloc( HYPRE_Bool, num_cols_offd );
 
 /*---------------------------------------------------------------------
 * If there exists no CommPkg for A, a CommPkg is generated using
 * equally load balanced partitionings
 *--------------------------------------------------------------------*/
 if( !comm_pkg )
 {
 hypre_MatvecCommPkgCreate( A );
 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
 }
 
 num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
 x_buf = hypre_CTAlloc( HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
 
 index = 0;
 for( i = 0; i < num_sends; i++ )
 {
 j = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
 for( ; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
 {
 x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j )];
 }
 }
 
 comm_handle = hypre_ParCSRCommHandleCreate_bool( 1, comm_pkg, x_buf, x_tmp );
 
 hypre_CSRMatrixBooleanMatvec( diag, alpha, x, beta, y );
 
 hypre_ParCSRCommHandleDestroy( comm_handle );
 
 if( num_cols_offd )
 {
 hypre_CSRMatrixBooleanMatvec( offd, alpha, x_tmp, 1, y );
 }
 
 hypre_TFree( x_buf );
 hypre_TFree( x_tmp );
 }
 
 HYPRE_Int hypre_CSRMatrixSum( hypre_CSRMatrix* A, HYPRE_Complex beta, hypre_CSRMatrix* B )
 {
 HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
 HYPRE_Int nrows_A = hypre_CSRMatrixNumRows( A );
 HYPRE_Int ncols_A = hypre_CSRMatrixNumCols( A );
 HYPRE_Complex* B_data = hypre_CSRMatrixData( B );
 HYPRE_Int* B_i = hypre_CSRMatrixI( B );
 HYPRE_Int* B_j = hypre_CSRMatrixJ( B );
 HYPRE_Int nrows_B = hypre_CSRMatrixNumRows( B );
 HYPRE_Int ncols_B = hypre_CSRMatrixNumCols( B );
 
 HYPRE_Int ia, j, pos;
 HYPRE_Int* marker;
 
 if( nrows_A != nrows_B || ncols_A != ncols_B )
 {
 return -1; /* error: incompatible matrix dimensions */
 }
 
 marker = hypre_CTAlloc( HYPRE_Int, ncols_A );
 for( ia = 0; ia < ncols_A; ia++ )
 {
 marker[ia] = -1;
 }
 
 for( ia = 0; ia < nrows_A; ia++ )
 {
 for( j = A_i[ia]; j < A_i[ia + 1]; j++ )
 {
 marker[A_j[j]] = j;
 }
 
 for( j = B_i[ia]; j < B_i[ia + 1]; j++ )
 {
 pos = marker[B_j[j]];
 if( pos < A_i[ia] )
 {
 return -2; /* error: found an entry in B that is not present in A */
 }
 A_data[pos] += beta * B_data[j];
 }
 }
 
 hypre_TFree( marker );
 return 0;
 }
 
 hypre_ParCSRMatrix* hypre_ParCSRMatrixAdd( hypre_ParCSRMatrix* A, hypre_ParCSRMatrix* B )
 {
 MPI_Comm comm = hypre_ParCSRMatrixComm( A );
 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
 hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
 HYPRE_Int* A_cmap = hypre_ParCSRMatrixColMapOffd( A );
 HYPRE_Int A_cmap_size = hypre_CSRMatrixNumCols( A_offd );
 hypre_CSRMatrix* B_diag = hypre_ParCSRMatrixDiag( B );
 hypre_CSRMatrix* B_offd = hypre_ParCSRMatrixOffd( B );
 HYPRE_Int* B_cmap = hypre_ParCSRMatrixColMapOffd( B );
 HYPRE_Int B_cmap_size = hypre_CSRMatrixNumCols( B_offd );
 hypre_ParCSRMatrix* C;
 hypre_CSRMatrix* C_diag;
 hypre_CSRMatrix* C_offd;
 HYPRE_Int* C_cmap;
 HYPRE_Int im;
 HYPRE_Int cmap_differ;
 
 /* Check if A_cmap and B_cmap are the same. */
 cmap_differ = 0;
 if( A_cmap_size != B_cmap_size )
 {
 cmap_differ = 1; /* A and B have different cmap_size */
 }
 else
 {
 for( im = 0; im < A_cmap_size; im++ )
 {
 if( A_cmap[im] != B_cmap[im] )
 {
 cmap_differ = 1; /* A and B have different cmap arrays */
 break;
 }
 }
 }
 
 if( cmap_differ == 0 )
 {
 /* A and B have the same column mapping for their off-diagonal blocks so
 we can sum the diagonal and off-diagonal blocks separately and reduce
 temporary memory usage. */
 
 /* Add diagonals, off-diagonals, copy cmap. */
 C_diag = hypre_CSRMatrixAdd( A_diag, B_diag );
 if( !C_diag )
 {
 return NULL; /* error: A_diag and B_diag have different dimensions */
 }
 C_offd = hypre_CSRMatrixAdd( A_offd, B_offd );
 if( !C_offd )
 {
 hypre_CSRMatrixDestroy( C_diag );
 return NULL; /* error: A_offd and B_offd have different dimensions */
 }
 /* copy A_cmap -> C_cmap */
 C_cmap = hypre_TAlloc( HYPRE_Int, A_cmap_size );
 for( im = 0; im < A_cmap_size; im++ )
 {
 C_cmap[im] = A_cmap[im];
 }
 
 C = hypre_ParCSRMatrixCreate( comm, hypre_ParCSRMatrixGlobalNumRows( A ),
 hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
 hypre_ParCSRMatrixColStarts( A ), hypre_CSRMatrixNumCols( C_offd ),
 hypre_CSRMatrixNumNonzeros( C_diag ), hypre_CSRMatrixNumNonzeros( C_offd ) );
 
 /* In C, destroy diag/offd (allocated by Create) and replace them with
 C_diag/C_offd. */
 hypre_CSRMatrixDestroy( hypre_ParCSRMatrixDiag( C ) );
 hypre_CSRMatrixDestroy( hypre_ParCSRMatrixOffd( C ) );
 hypre_ParCSRMatrixDiag( C ) = C_diag;
 hypre_ParCSRMatrixOffd( C ) = C_offd;
 
 hypre_ParCSRMatrixColMapOffd( C ) = C_cmap;
 }
 else
 {
 /* A and B have different column mappings for their off-diagonal blocks so
 we need to use the column maps to create full-width CSR matricies. */
 
 int ierr = 0;
 hypre_CSRMatrix* csr_A;
 hypre_CSRMatrix* csr_B;
 hypre_CSRMatrix* csr_C_temp;
 
 /* merge diag and off-diag portions of A */
 csr_A = hypre_MergeDiagAndOffd( A );
 
 /* merge diag and off-diag portions of B */
 csr_B = hypre_MergeDiagAndOffd( B );
 
 /* add A and B */
 csr_C_temp = hypre_CSRMatrixAdd( csr_A, csr_B );
 
 /* delete CSR versions of A and B */
 ierr += hypre_CSRMatrixDestroy( csr_A );
 ierr += hypre_CSRMatrixDestroy( csr_B );
 
 /* create a new empty ParCSR matrix to contain the sum */
 C = hypre_ParCSRMatrixCreate( hypre_ParCSRMatrixComm( A ), hypre_ParCSRMatrixGlobalNumRows( A ),
 hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
 hypre_ParCSRMatrixColStarts( A ), 0, 0, 0 );
 
 /* split C into diag and off-diag portions */
 /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the
 number of columns in csr_C_temp which is the global number of columns
 in A and B. This does not scale well. */
 ierr += GenerateDiagAndOffd( csr_C_temp, C, hypre_ParCSRMatrixFirstColDiag( A ),
 hypre_ParCSRMatrixLastColDiag( A ) );
 
 /* delete CSR version of C */
 ierr += hypre_CSRMatrixDestroy( csr_C_temp );
 }
 
 /* hypre_ParCSRMatrixSetNumNonzeros(A); */
 
 /* Make sure that the first entry in each row is the diagonal one. */
 hypre_CSRMatrixReorder( hypre_ParCSRMatrixDiag( C ) );
 
 /* C owns diag, offd, and cmap. */
 hypre_ParCSRMatrixSetDataOwner( C, 1 );
 /* C does not own row and column starts. */
 hypre_ParCSRMatrixSetRowStartsOwner( C, 0 );
 hypre_ParCSRMatrixSetColStartsOwner( C, 0 );
 
 return C;
 }
 
 HYPRE_Int hypre_ParCSRMatrixSum( hypre_ParCSRMatrix* A, HYPRE_Complex beta, hypre_ParCSRMatrix* B )
 {
 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
 hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
 hypre_CSRMatrix* B_diag = hypre_ParCSRMatrixDiag( B );
 hypre_CSRMatrix* B_offd = hypre_ParCSRMatrixOffd( B );
 HYPRE_Int error;
 
 error = hypre_CSRMatrixSum( A_diag, beta, B_diag );
 error = error ? error : hypre_CSRMatrixSum( A_offd, beta, B_offd );
 
 return error;
 }
 
 HYPRE_Int hypre_CSRMatrixSetConstantValues( hypre_CSRMatrix* A, HYPRE_Complex value )
 {
 HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
 HYPRE_Int A_nnz = hypre_CSRMatrixNumNonzeros( A );
 HYPRE_Int ia;
 
 for( ia = 0; ia < A_nnz; ia++ )
 {
 A_data[ia] = value;
 }
 
 return 0;
 }
 
 HYPRE_Int hypre_ParCSRMatrixSetConstantValues( hypre_ParCSRMatrix* A, HYPRE_Complex value )
 {
 hypre_CSRMatrixSetConstantValues( hypre_ParCSRMatrixDiag( A ), value );
 hypre_CSRMatrixSetConstantValues( hypre_ParCSRMatrixOffd( A ), value );
 
 return 0;
 }
 
} // namespace internal
 
} // namespace moab
 
#endif // MOAB_HAVE_MPI