Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
hypre_parcsr.cpp
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 #include "hypre_parcsr.hpp"
00013 
00014 #ifdef MOAB_HAVE_MPI
00015 
00016 #include <limits>
00017 
00018 namespace moab
00019 {
00020 
00021 namespace internal
00022 {
00023 
00024     /*--------------------------------------------------------------------------
00025      *                        A*X = B style elimination
00026      *--------------------------------------------------------------------------*/
00027 
00028     /*
00029         Function:  hypre_CSRMatrixEliminateAXB
00030 
00031         Eliminate the rows and columns of A corresponding to the
00032         given sorted (!) list of rows. Put I on the eliminated diagonal,
00033         subtract columns times X from B.
00034     */
00035     void hypre_CSRMatrixEliminateAXB( hypre_CSRMatrix* A,
00036                                       HYPRE_Int nrows_to_eliminate,
00037                                       HYPRE_Int* rows_to_eliminate,
00038                                       hypre_Vector* X,
00039                                       hypre_Vector* B )
00040     {
00041         HYPRE_Int i, j;
00042         HYPRE_Int irow, jcol, ibeg, iend, pos;
00043         HYPRE_Real a;
00044 
00045         HYPRE_Int* Ai     = hypre_CSRMatrixI( A );
00046         HYPRE_Int* Aj     = hypre_CSRMatrixJ( A );
00047         HYPRE_Real* Adata = hypre_CSRMatrixData( A );
00048         HYPRE_Int nrows   = hypre_CSRMatrixNumRows( A );
00049 
00050         HYPRE_Real* Xdata = hypre_VectorData( X );
00051         HYPRE_Real* Bdata = hypre_VectorData( B );
00052 
00053         /* eliminate the columns */
00054         for( i = 0; i < nrows; i++ )
00055         {
00056             ibeg = Ai[i];
00057             iend = Ai[i + 1];
00058             for( j = ibeg; j < iend; j++ )
00059             {
00060                 jcol = Aj[j];
00061                 pos  = hypre_BinarySearch( rows_to_eliminate, jcol, nrows_to_eliminate );
00062                 if( pos != -1 )
00063                 {
00064                     a        = Adata[j];
00065                     Adata[j] = 0.0;
00066                     Bdata[i] -= a * Xdata[jcol];
00067                 }
00068             }
00069         }
00070 
00071         /* remove the rows and set the diagonal equal to 1 */
00072         for( i = 0; i < nrows_to_eliminate; i++ )
00073         {
00074             irow = rows_to_eliminate[i];
00075             ibeg = Ai[irow];
00076             iend = Ai[irow + 1];
00077             for( j = ibeg; j < iend; j++ )
00078             {
00079                 if( Aj[j] == irow )
00080                 {
00081                     Adata[j] = 1.0;
00082                 }
00083                 else
00084                 {
00085                     Adata[j] = 0.0;
00086                 }
00087             }
00088         }
00089     }
00090 
00091     /*
00092         Function:  hypre_CSRMatrixEliminateOffdColsAXB
00093 
00094         Eliminate the given sorted (!) list of columns of A, subtract them from B.
00095     */
00096     void hypre_CSRMatrixEliminateOffdColsAXB( hypre_CSRMatrix* A,
00097                                               HYPRE_Int ncols_to_eliminate,
00098                                               HYPRE_Int* eliminate_cols,
00099                                               HYPRE_Real* eliminate_coefs,
00100                                               hypre_Vector* B )
00101     {
00102         HYPRE_Int i, j;
00103         HYPRE_Int ibeg, iend, pos;
00104         HYPRE_Real a;
00105 
00106         HYPRE_Int* Ai     = hypre_CSRMatrixI( A );
00107         HYPRE_Int* Aj     = hypre_CSRMatrixJ( A );
00108         HYPRE_Real* Adata = hypre_CSRMatrixData( A );
00109         HYPRE_Int nrows   = hypre_CSRMatrixNumRows( A );
00110 
00111         HYPRE_Real* Bdata = hypre_VectorData( B );
00112 
00113         for( i = 0; i < nrows; i++ )
00114         {
00115             ibeg = Ai[i];
00116             iend = Ai[i + 1];
00117             for( j = ibeg; j < iend; j++ )
00118             {
00119                 pos = hypre_BinarySearch( eliminate_cols, Aj[j], ncols_to_eliminate );
00120                 if( pos != -1 )
00121                 {
00122                     a        = Adata[j];
00123                     Adata[j] = 0.0;
00124                     Bdata[i] -= a * eliminate_coefs[pos];
00125                 }
00126             }
00127         }
00128     }
00129 
00130     /*
00131         Function:  hypre_CSRMatrixEliminateOffdRowsAXB
00132 
00133         Eliminate (zero) the given list of rows of A.
00134     */
00135     void hypre_CSRMatrixEliminateOffdRowsAXB( hypre_CSRMatrix* A,
00136                                               HYPRE_Int nrows_to_eliminate,
00137                                               HYPRE_Int* rows_to_eliminate )
00138     {
00139         HYPRE_Int* Ai     = hypre_CSRMatrixI( A );
00140         HYPRE_Real* Adata = hypre_CSRMatrixData( A );
00141 
00142         HYPRE_Int i, j;
00143         HYPRE_Int irow, ibeg, iend;
00144 
00145         for( i = 0; i < nrows_to_eliminate; i++ )
00146         {
00147             irow = rows_to_eliminate[i];
00148             ibeg = Ai[irow];
00149             iend = Ai[irow + 1];
00150             for( j = ibeg; j < iend; j++ )
00151             {
00152                 Adata[j] = 0.0;
00153             }
00154         }
00155     }
00156 
00157     /*
00158         Function:  hypre_ParCSRMatrixEliminateAXB
00159 
00160         This function eliminates the global rows and columns of a matrix
00161         A corresponding to given lists of sorted (!) local row numbers,
00162         so that the solution to the system A*X = B is X_b for the given rows.
00163 
00164         The elimination is done as follows:
00165 
00166                                             (input)                  (output)
00167 
00168                                     / A_ii | A_ib \          / A_ii |  0   \
00169                             A = | -----+----- |   --->   | -----+----- |
00170                                     \ A_bi | A_bb /          \   0  |  I   /
00171 
00172                                                     / X_i \          / X_i \
00173                                             X = | --- |   --->   | --- |  (no change)
00174                                                     \ X_b /          \ X_b /
00175 
00176                                                     / B_i \          / B_i - A_ib * X_b \
00177                                             B = | --- |   --->   | ---------------- |
00178                                                     \ B_b /          \        X_b       /
00179 
00180     */
00181     void hypre_ParCSRMatrixEliminateAXB( hypre_ParCSRMatrix* A,
00182                                          HYPRE_Int num_rowscols_to_elim,
00183                                          HYPRE_Int* rowscols_to_elim,
00184                                          hypre_ParVector* X,
00185                                          hypre_ParVector* B )
00186     {
00187         hypre_CSRMatrix* diag = hypre_ParCSRMatrixDiag( A );
00188         hypre_CSRMatrix* offd = hypre_ParCSRMatrixOffd( A );
00189         HYPRE_Int diag_nrows  = hypre_CSRMatrixNumRows( diag );
00190         HYPRE_Int offd_ncols  = hypre_CSRMatrixNumCols( offd );
00191 
00192         hypre_Vector* Xlocal = hypre_ParVectorLocalVector( X );
00193         hypre_Vector* Blocal = hypre_ParVectorLocalVector( B );
00194 
00195         HYPRE_Real* Bdata = hypre_VectorData( Blocal );
00196         HYPRE_Real* Xdata = hypre_VectorData( Xlocal );
00197 
00198         HYPRE_Int num_offd_cols_to_elim;
00199         HYPRE_Int* offd_cols_to_elim;
00200         HYPRE_Real* eliminate_coefs;
00201 
00202         /* figure out which offd cols should be eliminated and with what coef */
00203         hypre_ParCSRCommHandle* comm_handle;
00204         hypre_ParCSRCommPkg* comm_pkg;
00205         HYPRE_Int num_sends;
00206         HYPRE_Int index, start;
00207         HYPRE_Int i, j, k, irow;
00208 
00209         HYPRE_Real* eliminate_row = hypre_CTAlloc( HYPRE_Real, diag_nrows );
00210         HYPRE_Real* eliminate_col = hypre_CTAlloc( HYPRE_Real, offd_ncols );
00211         HYPRE_Real *buf_data, coef;
00212 
00213         /* make sure A has a communication package */
00214         comm_pkg = hypre_ParCSRMatrixCommPkg( A );
00215         if( !comm_pkg )
00216         {
00217             hypre_MatvecCommPkgCreate( A );
00218             comm_pkg = hypre_ParCSRMatrixCommPkg( A );
00219         }
00220 
00221         /* HACK: rows that shouldn't be eliminated are marked with quiet NaN;
00222                 those that should are set to the boundary value from X; this is to
00223                 avoid sending complex type (int+double) or communicating twice. */
00224         for( i = 0; i < diag_nrows; i++ )
00225         {
00226             eliminate_row[i] = std::numeric_limits< HYPRE_Real >::quiet_NaN();
00227         }
00228         for( i = 0; i < num_rowscols_to_elim; i++ )
00229         {
00230             irow                = rowscols_to_elim[i];
00231             eliminate_row[irow] = Xdata[irow];
00232         }
00233 
00234         /* use a Matvec communication pattern to find (in eliminate_col)
00235                 which of the local offd columns are to be eliminated */
00236         num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
00237         buf_data  = hypre_CTAlloc( HYPRE_Real, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
00238         index     = 0;
00239         for( i = 0; i < num_sends; i++ )
00240         {
00241             start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
00242             for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
00243             {
00244                 k                 = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
00245                 buf_data[index++] = eliminate_row[k];
00246             }
00247         }
00248         comm_handle = hypre_ParCSRCommHandleCreate( 1, comm_pkg, buf_data, eliminate_col );
00249 
00250         /* do sequential part of the elimination while stuff is getting sent */
00251         hypre_CSRMatrixEliminateAXB( diag, num_rowscols_to_elim, rowscols_to_elim, Xlocal, Blocal );
00252 
00253         /* finish the communication */
00254         hypre_ParCSRCommHandleDestroy( comm_handle );
00255 
00256         /* received eliminate_col[], count offd columns to eliminate */
00257         num_offd_cols_to_elim = 0;
00258         for( i = 0; i < offd_ncols; i++ )
00259         {
00260             coef = eliminate_col[i];
00261             if( coef == coef )  // test for NaN
00262             {
00263                 num_offd_cols_to_elim++;
00264             }
00265         }
00266 
00267         offd_cols_to_elim = hypre_CTAlloc( HYPRE_Int, num_offd_cols_to_elim );
00268         eliminate_coefs   = hypre_CTAlloc( HYPRE_Real, num_offd_cols_to_elim );
00269 
00270         /* get a list of offd column indices and coefs */
00271         num_offd_cols_to_elim = 0;
00272         for( i = 0; i < offd_ncols; i++ )
00273         {
00274             coef = eliminate_col[i];
00275             if( coef == coef )  // test for NaN
00276             {
00277                 offd_cols_to_elim[num_offd_cols_to_elim] = i;
00278                 eliminate_coefs[num_offd_cols_to_elim]   = coef;
00279                 num_offd_cols_to_elim++;
00280             }
00281         }
00282 
00283         hypre_TFree( buf_data );
00284         hypre_TFree( eliminate_row );
00285         hypre_TFree( eliminate_col );
00286 
00287         /* eliminate the off-diagonal part */
00288         hypre_CSRMatrixEliminateOffdColsAXB( offd, num_offd_cols_to_elim, offd_cols_to_elim, eliminate_coefs, Blocal );
00289 
00290         hypre_CSRMatrixEliminateOffdRowsAXB( offd, num_rowscols_to_elim, rowscols_to_elim );
00291 
00292         /* set boundary values in the rhs */
00293         for( int m = 0; m < num_rowscols_to_elim; m++ )
00294         {
00295             irow        = rowscols_to_elim[m];
00296             Bdata[irow] = Xdata[irow];
00297         }
00298 
00299         hypre_TFree( offd_cols_to_elim );
00300         hypre_TFree( eliminate_coefs );
00301     }
00302 
00303     /*--------------------------------------------------------------------------
00304      *                        (A + Ae) style elimination
00305      *--------------------------------------------------------------------------*/
00306 
00307     /*
00308         Function:  hypre_CSRMatrixElimCreate
00309 
00310         Prepare the Ae matrix: count nnz, initialize I, allocate J and data.
00311     */
00312     void hypre_CSRMatrixElimCreate( hypre_CSRMatrix* A,
00313                                     hypre_CSRMatrix* Ae,
00314                                     HYPRE_Int nrows,
00315                                     HYPRE_Int* rows,
00316                                     HYPRE_Int ncols,
00317                                     HYPRE_Int* cols,
00318                                     HYPRE_Int* col_mark )
00319     {
00320         HYPRE_Int i, j, col;
00321         HYPRE_Int A_beg, A_end;
00322 
00323         HYPRE_Int* A_i   = hypre_CSRMatrixI( A );
00324         HYPRE_Int* A_j   = hypre_CSRMatrixJ( A );
00325         HYPRE_Int A_rows = hypre_CSRMatrixNumRows( A );
00326 
00327         hypre_CSRMatrixI( Ae ) = hypre_TAlloc( HYPRE_Int, A_rows + 1 );
00328 
00329         HYPRE_Int* Ae_i = hypre_CSRMatrixI( Ae );
00330         HYPRE_Int nnz   = 0;
00331 
00332         for( i = 0; i < A_rows; i++ )
00333         {
00334             Ae_i[i] = nnz;
00335 
00336             A_beg = A_i[i];
00337             A_end = A_i[i + 1];
00338 
00339             if( hypre_BinarySearch( rows, i, nrows ) >= 0 )
00340             {
00341                 /* full row */
00342                 nnz += A_end - A_beg;
00343 
00344                 if( col_mark )
00345                 {
00346                     for( j = A_beg; j < A_end; j++ )
00347                     {
00348                         col_mark[A_j[j]] = 1;
00349                     }
00350                 }
00351             }
00352             else
00353             {
00354                 /* count columns */
00355                 for( j = A_beg; j < A_end; j++ )
00356                 {
00357                     col = A_j[j];
00358                     if( hypre_BinarySearch( cols, col, ncols ) >= 0 )
00359                     {
00360                         nnz++;
00361                         if( col_mark )
00362                         {
00363                             col_mark[col] = 1;
00364                         }
00365                     }
00366                 }
00367             }
00368         }
00369         Ae_i[A_rows] = nnz;
00370 
00371         hypre_CSRMatrixJ( Ae )           = hypre_TAlloc( HYPRE_Int, nnz );
00372         hypre_CSRMatrixData( Ae )        = hypre_TAlloc( HYPRE_Real, nnz );
00373         hypre_CSRMatrixNumNonzeros( Ae ) = nnz;
00374     }
00375 
00376     /*
00377         Function:  hypre_CSRMatrixEliminateRowsCols
00378 
00379         Eliminate rows and columns of A, store eliminated values in Ae.
00380         If 'diag' is nonzero, the eliminated diagonal of A is set to identity.
00381         If 'col_remap' is not NULL it specifies renumbering of columns of Ae.
00382     */
00383     void hypre_CSRMatrixEliminateRowsCols( hypre_CSRMatrix* A,
00384                                            hypre_CSRMatrix* Ae,
00385                                            HYPRE_Int nrows,
00386                                            HYPRE_Int* rows,
00387                                            HYPRE_Int ncols,
00388                                            HYPRE_Int* cols,
00389                                            int diag,
00390                                            HYPRE_Int* col_remap )
00391     {
00392         HYPRE_Int i, j, k, col;
00393         HYPRE_Int A_beg, Ae_beg, A_end;
00394         HYPRE_Real a;
00395 
00396         HYPRE_Int* A_i     = hypre_CSRMatrixI( A );
00397         HYPRE_Int* A_j     = hypre_CSRMatrixJ( A );
00398         HYPRE_Real* A_data = hypre_CSRMatrixData( A );
00399         HYPRE_Int A_rows   = hypre_CSRMatrixNumRows( A );
00400 
00401         HYPRE_Int* Ae_i     = hypre_CSRMatrixI( Ae );
00402         HYPRE_Int* Ae_j     = hypre_CSRMatrixJ( Ae );
00403         HYPRE_Real* Ae_data = hypre_CSRMatrixData( Ae );
00404 
00405         for( i = 0; i < A_rows; i++ )
00406         {
00407             A_beg  = A_i[i];
00408             A_end  = A_i[i + 1];
00409             Ae_beg = Ae_i[i];
00410 
00411             if( hypre_BinarySearch( rows, i, nrows ) >= 0 )
00412             {
00413                 /* eliminate row */
00414                 for( j = A_beg, k = Ae_beg; j < A_end; j++, k++ )
00415                 {
00416                     col        = A_j[j];
00417                     Ae_j[k]    = col_remap ? col_remap[col] : col;
00418                     a          = ( diag && col == i ) ? 1.0 : 0.0;
00419                     Ae_data[k] = A_data[j] - a;
00420                     A_data[j]  = a;
00421                 }
00422             }
00423             else
00424             {
00425                 /* eliminate columns */
00426                 for( j = A_beg, k = Ae_beg; j < A_end; j++ )
00427                 {
00428                     col = A_j[j];
00429                     if( hypre_BinarySearch( cols, col, ncols ) >= 0 )
00430                     {
00431                         Ae_j[k]    = col_remap ? col_remap[col] : col;
00432                         Ae_data[k] = A_data[j];
00433                         A_data[j]  = 0.0;
00434                         k++;
00435                     }
00436                 }
00437             }
00438         }
00439     }
00440 
00441     /*
00442         Function:  hypre_ParCSRMatrixEliminateAAe
00443 
00444                                             (input)                  (output)
00445 
00446                                     / A_ii | A_ib \          / A_ii |  0   \
00447                             A = | -----+----- |   --->   | -----+----- |
00448                                     \ A_bi | A_bb /          \   0  |  I   /
00449 
00450 
00451                                                                                     /   0  |   A_ib
00452        \ Ae = | -----+--------- | \ A_bi | A_bb - I /
00453 
00454     */
00455     void hypre_ParCSRMatrixEliminateAAe( hypre_ParCSRMatrix* A,
00456                                          hypre_ParCSRMatrix** Ae,
00457                                          HYPRE_Int num_rowscols_to_elim,
00458                                          HYPRE_Int* rowscols_to_elim )
00459     {
00460         HYPRE_Int i, j, k;
00461 
00462         hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
00463         hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
00464         HYPRE_Int A_diag_nrows  = hypre_CSRMatrixNumRows( A_diag );
00465         HYPRE_Int A_offd_ncols  = hypre_CSRMatrixNumCols( A_offd );
00466 
00467         *Ae = hypre_ParCSRMatrixCreate( hypre_ParCSRMatrixComm( A ), hypre_ParCSRMatrixGlobalNumRows( A ),
00468                                         hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
00469                                         hypre_ParCSRMatrixColStarts( A ), 0, 0, 0 );
00470 
00471         hypre_ParCSRMatrixSetRowStartsOwner( *Ae, 0 );
00472         hypre_ParCSRMatrixSetColStartsOwner( *Ae, 0 );
00473 
00474         hypre_CSRMatrix* Ae_diag = hypre_ParCSRMatrixDiag( *Ae );
00475         hypre_CSRMatrix* Ae_offd = hypre_ParCSRMatrixOffd( *Ae );
00476         HYPRE_Int Ae_offd_ncols;
00477 
00478         HYPRE_Int num_offd_cols_to_elim;
00479         HYPRE_Int* offd_cols_to_elim;
00480 
00481         HYPRE_Int* A_col_map_offd = hypre_ParCSRMatrixColMapOffd( A );
00482         HYPRE_Int* Ae_col_map_offd;
00483 
00484         HYPRE_Int* col_mark;
00485         HYPRE_Int* col_remap;
00486 
00487         /* figure out which offd cols should be eliminated */
00488         {
00489             hypre_ParCSRCommHandle* comm_handle;
00490             hypre_ParCSRCommPkg* comm_pkg;
00491             HYPRE_Int num_sends, *int_buf_data;
00492             HYPRE_Int index, start;
00493 
00494             HYPRE_Int* eliminate_row = hypre_CTAlloc( HYPRE_Int, A_diag_nrows );
00495             HYPRE_Int* eliminate_col = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
00496 
00497             /* make sure A has a communication package */
00498             comm_pkg = hypre_ParCSRMatrixCommPkg( A );
00499             if( !comm_pkg )
00500             {
00501                 hypre_MatvecCommPkgCreate( A );
00502                 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
00503             }
00504 
00505             /* which of the local rows are to be eliminated */
00506             for( i = 0; i < A_diag_nrows; i++ )
00507             {
00508                 eliminate_row[i] = 0;
00509             }
00510             for( i = 0; i < num_rowscols_to_elim; i++ )
00511             {
00512                 eliminate_row[rowscols_to_elim[i]] = 1;
00513             }
00514 
00515             /* use a Matvec communication pattern to find (in eliminate_col)
00516                 which of the local offd columns are to be eliminated */
00517             num_sends    = hypre_ParCSRCommPkgNumSends( comm_pkg );
00518             int_buf_data = hypre_CTAlloc( HYPRE_Int, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
00519             index        = 0;
00520             for( i = 0; i < num_sends; i++ )
00521             {
00522                 start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
00523                 for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
00524                 {
00525                     k                     = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
00526                     int_buf_data[index++] = eliminate_row[k];
00527                 }
00528             }
00529             comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data, eliminate_col );
00530 
00531             /* eliminate diagonal part, overlapping it with communication */
00532             hypre_CSRMatrixElimCreate( A_diag, Ae_diag, num_rowscols_to_elim, rowscols_to_elim, num_rowscols_to_elim,
00533                                        rowscols_to_elim, NULL );
00534 
00535             hypre_CSRMatrixEliminateRowsCols( A_diag, Ae_diag, num_rowscols_to_elim, rowscols_to_elim,
00536                                               num_rowscols_to_elim, rowscols_to_elim, 1, NULL );
00537             hypre_CSRMatrixReorder( Ae_diag );
00538 
00539             /* finish the communication */
00540             hypre_ParCSRCommHandleDestroy( comm_handle );
00541 
00542             /* received eliminate_col[], count offd columns to eliminate */
00543             num_offd_cols_to_elim = 0;
00544             for( i = 0; i < A_offd_ncols; i++ )
00545             {
00546                 if( eliminate_col[i] )
00547                 {
00548                     num_offd_cols_to_elim++;
00549                 }
00550             }
00551 
00552             offd_cols_to_elim = hypre_CTAlloc( HYPRE_Int, num_offd_cols_to_elim );
00553 
00554             /* get a list of offd column indices and coefs */
00555             num_offd_cols_to_elim = 0;
00556             for( i = 0; i < A_offd_ncols; i++ )
00557             {
00558                 if( eliminate_col[i] )
00559                 {
00560                     offd_cols_to_elim[num_offd_cols_to_elim++] = i;
00561                 }
00562             }
00563 
00564             hypre_TFree( int_buf_data );
00565             hypre_TFree( eliminate_row );
00566             hypre_TFree( eliminate_col );
00567         }
00568 
00569         /* eliminate the off-diagonal part */
00570         col_mark  = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
00571         col_remap = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
00572 
00573         hypre_CSRMatrixElimCreate( A_offd, Ae_offd, num_rowscols_to_elim, rowscols_to_elim, num_offd_cols_to_elim,
00574                                    offd_cols_to_elim, col_mark );
00575 
00576         for( i = k = 0; i < A_offd_ncols; i++ )
00577         {
00578             if( col_mark[i] )
00579             {
00580                 col_remap[i] = k++;
00581             }
00582         }
00583 
00584         hypre_CSRMatrixEliminateRowsCols( A_offd, Ae_offd, num_rowscols_to_elim, rowscols_to_elim,
00585                                           num_offd_cols_to_elim, offd_cols_to_elim, 0, col_remap );
00586 
00587         /* create col_map_offd for Ae */
00588         Ae_offd_ncols = 0;
00589         for( i = 0; i < A_offd_ncols; i++ )
00590         {
00591             if( col_mark[i] )
00592             {
00593                 Ae_offd_ncols++;
00594             }
00595         }
00596 
00597         Ae_col_map_offd = hypre_CTAlloc( HYPRE_Int, Ae_offd_ncols );
00598 
00599         Ae_offd_ncols = 0;
00600         for( i = 0; i < A_offd_ncols; i++ )
00601         {
00602             if( col_mark[i] )
00603             {
00604                 Ae_col_map_offd[Ae_offd_ncols++] = A_col_map_offd[i];
00605             }
00606         }
00607 
00608         hypre_ParCSRMatrixColMapOffd( *Ae ) = Ae_col_map_offd;
00609         hypre_CSRMatrixNumCols( Ae_offd )   = Ae_offd_ncols;
00610 
00611         hypre_TFree( col_remap );
00612         hypre_TFree( col_mark );
00613         hypre_TFree( offd_cols_to_elim );
00614 
00615         hypre_ParCSRMatrixSetNumNonzeros( *Ae );
00616         hypre_MatvecCommPkgCreate( *Ae );
00617     }
00618 
00619     /*--------------------------------------------------------------------------
00620      *                               Split
00621      *--------------------------------------------------------------------------*/
00622 
00623     void hypre_CSRMatrixSplit( hypre_CSRMatrix* A,
00624                                HYPRE_Int nr,
00625                                HYPRE_Int nc,
00626                                HYPRE_Int* row_block_num,
00627                                HYPRE_Int* col_block_num,
00628                                hypre_CSRMatrix** blocks )
00629     {
00630         HYPRE_Int i, j, k, bi, bj;
00631 
00632         HYPRE_Int* A_i        = hypre_CSRMatrixI( A );
00633         HYPRE_Int* A_j        = hypre_CSRMatrixJ( A );
00634         HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
00635 
00636         HYPRE_Int A_rows = hypre_CSRMatrixNumRows( A );
00637         HYPRE_Int A_cols = hypre_CSRMatrixNumCols( A );
00638 
00639         HYPRE_Int* num_rows = hypre_CTAlloc( HYPRE_Int, nr );
00640         HYPRE_Int* num_cols = hypre_CTAlloc( HYPRE_Int, nc );
00641 
00642         HYPRE_Int* block_row = hypre_TAlloc( HYPRE_Int, A_rows );
00643         HYPRE_Int* block_col = hypre_TAlloc( HYPRE_Int, A_cols );
00644 
00645         for( i = 0; i < A_rows; i++ )
00646         {
00647             block_row[i] = num_rows[row_block_num[i]]++;
00648         }
00649         for( j = 0; j < A_cols; j++ )
00650         {
00651             block_col[j] = num_cols[col_block_num[j]]++;
00652         }
00653 
00654         /* allocate the blocks */
00655         for( i = 0; i < nr; i++ )
00656         {
00657             for( j = 0; j < nc; j++ )
00658             {
00659                 hypre_CSRMatrix* B    = hypre_CSRMatrixCreate( num_rows[i], num_cols[j], 0 );
00660                 hypre_CSRMatrixI( B ) = hypre_CTAlloc( HYPRE_Int, num_rows[i] + 1 );
00661                 blocks[i * nc + j]    = B;
00662             }
00663         }
00664 
00665         /* count block row nnz */
00666         for( i = 0; i < A_rows; i++ )
00667         {
00668             bi = row_block_num[i];
00669             for( j = A_i[i]; j < A_i[i + 1]; j++ )
00670             {
00671                 bj                 = col_block_num[A_j[j]];
00672                 hypre_CSRMatrix* B = blocks[bi * nc + bj];
00673                 hypre_CSRMatrixI( B )[block_row[i] + 1]++;
00674             }
00675         }
00676 
00677         /* count block nnz */
00678         for( k = 0; k < nr * nc; k++ )
00679         {
00680             hypre_CSRMatrix* B = blocks[k];
00681             HYPRE_Int* B_i     = hypre_CSRMatrixI( B );
00682 
00683             HYPRE_Int nnz = 0, rs;
00684             for( int m = 1; m <= hypre_CSRMatrixNumRows( B ); m++ )
00685             {
00686                 rs = B_i[m], B_i[m] = nnz, nnz += rs;
00687             }
00688 
00689             hypre_CSRMatrixJ( B )           = hypre_TAlloc( HYPRE_Int, nnz );
00690             hypre_CSRMatrixData( B )        = hypre_TAlloc( HYPRE_Complex, nnz );
00691             hypre_CSRMatrixNumNonzeros( B ) = nnz;
00692         }
00693 
00694         /* populate blocks */
00695         for( i = 0; i < A_rows; i++ )
00696         {
00697             bi = row_block_num[i];
00698             for( j = A_i[i]; j < A_i[i + 1]; j++ )
00699             {
00700                 k                              = A_j[j];
00701                 bj                             = col_block_num[k];
00702                 hypre_CSRMatrix* B             = blocks[bi * nc + bj];
00703                 HYPRE_Int* bii                 = hypre_CSRMatrixI( B ) + block_row[i] + 1;
00704                 hypre_CSRMatrixJ( B )[*bii]    = block_col[k];
00705                 hypre_CSRMatrixData( B )[*bii] = A_data[j];
00706                 ( *bii )++;
00707             }
00708         }
00709 
00710         hypre_TFree( block_col );
00711         hypre_TFree( block_row );
00712 
00713         hypre_TFree( num_cols );
00714         hypre_TFree( num_rows );
00715     }
00716 
00717     void hypre_ParCSRMatrixSplit( hypre_ParCSRMatrix* A,
00718                                   HYPRE_Int nr,
00719                                   HYPRE_Int nc,
00720                                   hypre_ParCSRMatrix** blocks,
00721                                   int interleaved_rows,
00722                                   int interleaved_cols )
00723     {
00724         HYPRE_Int i, j, k;
00725 
00726         MPI_Comm comm = hypre_ParCSRMatrixComm( A );
00727 
00728         hypre_CSRMatrix* Adiag = hypre_ParCSRMatrixDiag( A );
00729         hypre_CSRMatrix* Aoffd = hypre_ParCSRMatrixOffd( A );
00730 
00731         HYPRE_Int global_rows = hypre_ParCSRMatrixGlobalNumRows( A );
00732         HYPRE_Int global_cols = hypre_ParCSRMatrixGlobalNumCols( A );
00733 
00734         HYPRE_Int local_rows = hypre_CSRMatrixNumRows( Adiag );
00735         HYPRE_Int local_cols = hypre_CSRMatrixNumCols( Adiag );
00736         HYPRE_Int offd_cols  = hypre_CSRMatrixNumCols( Aoffd );
00737 
00738         hypre_assert( local_rows % nr == 0 && local_cols % nc == 0 );
00739         hypre_assert( global_rows % nr == 0 && global_cols % nc == 0 );
00740 
00741         HYPRE_Int block_rows = local_rows / nr;
00742         HYPRE_Int block_cols = local_cols / nc;
00743         HYPRE_Int num_blocks = nr * nc;
00744 
00745         /* mark local rows and columns with block number */
00746         HYPRE_Int* row_block_num = hypre_TAlloc( HYPRE_Int, local_rows );
00747         HYPRE_Int* col_block_num = hypre_TAlloc( HYPRE_Int, local_cols );
00748 
00749         for( i = 0; i < local_rows; i++ )
00750         {
00751             row_block_num[i] = interleaved_rows ? ( i % nr ) : ( i / block_rows );
00752         }
00753         for( i = 0; i < local_cols; i++ )
00754         {
00755             col_block_num[i] = interleaved_cols ? ( i % nc ) : ( i / block_cols );
00756         }
00757 
00758         /* determine the block numbers for offd columns */
00759         HYPRE_Int* offd_col_block_num = hypre_TAlloc( HYPRE_Int, offd_cols );
00760         hypre_ParCSRCommHandle* comm_handle;
00761         HYPRE_Int* int_buf_data;
00762         {
00763             /* make sure A has a communication package */
00764             hypre_ParCSRCommPkg* comm_pkg = hypre_ParCSRMatrixCommPkg( A );
00765             if( !comm_pkg )
00766             {
00767                 hypre_MatvecCommPkgCreate( A );
00768                 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
00769             }
00770 
00771             /* calculate the final global column numbers for each block */
00772             HYPRE_Int* count            = hypre_CTAlloc( HYPRE_Int, nc );
00773             HYPRE_Int* block_global_col = hypre_TAlloc( HYPRE_Int, local_cols );
00774             HYPRE_Int first_col         = hypre_ParCSRMatrixFirstColDiag( A ) / nc;
00775             for( i = 0; i < local_cols; i++ )
00776             {
00777                 block_global_col[i] = first_col + count[col_block_num[i]]++;
00778             }
00779             hypre_TFree( count );
00780 
00781             /* use a Matvec communication pattern to determine offd_col_block_num */
00782             HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
00783             int_buf_data        = hypre_CTAlloc( HYPRE_Int, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
00784             HYPRE_Int start, index = 0;
00785             for( i = 0; i < num_sends; i++ )
00786             {
00787                 start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
00788                 for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
00789                 {
00790                     k                     = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
00791                     int_buf_data[index++] = col_block_num[k] + nc * block_global_col[k];
00792                 }
00793             }
00794             hypre_TFree( block_global_col );
00795 
00796             comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data, offd_col_block_num );
00797         }
00798 
00799         /* create the block matrices */
00800         HYPRE_Int num_procs = 1;
00801         if( !hypre_ParCSRMatrixAssumedPartition( A ) )
00802         {
00803             hypre_MPI_Comm_size( comm, &num_procs );
00804         }
00805 
00806         HYPRE_Int* row_starts = hypre_TAlloc( HYPRE_Int, num_procs + 1 );
00807         HYPRE_Int* col_starts = hypre_TAlloc( HYPRE_Int, num_procs + 1 );
00808         for( i = 0; i <= num_procs; i++ )
00809         {
00810             row_starts[i] = hypre_ParCSRMatrixRowStarts( A )[i] / nr;
00811             col_starts[i] = hypre_ParCSRMatrixColStarts( A )[i] / nc;
00812         }
00813 
00814         for( i = 0; i < num_blocks; i++ )
00815         {
00816             blocks[i] =
00817                 hypre_ParCSRMatrixCreate( comm, global_rows / nr, global_cols / nc, row_starts, col_starts, 0, 0, 0 );
00818         }
00819 
00820         /* split diag part */
00821         hypre_CSRMatrix** csr_blocks = hypre_TAlloc( hypre_CSRMatrix*, nr * nc );
00822         hypre_CSRMatrixSplit( Adiag, nr, nc, row_block_num, col_block_num, csr_blocks );
00823 
00824         for( i = 0; i < num_blocks; i++ )
00825         {
00826             hypre_TFree( hypre_ParCSRMatrixDiag( blocks[i] ) );
00827             hypre_ParCSRMatrixDiag( blocks[i] ) = csr_blocks[i];
00828         }
00829 
00830         /* finish communication, receive offd_col_block_num */
00831         hypre_ParCSRCommHandleDestroy( comm_handle );
00832         hypre_TFree( int_buf_data );
00833 
00834         /* decode global offd column numbers */
00835         HYPRE_Int* offd_global_col = hypre_TAlloc( HYPRE_Int, offd_cols );
00836         for( i = 0; i < offd_cols; i++ )
00837         {
00838             offd_global_col[i] = offd_col_block_num[i] / nc;
00839             offd_col_block_num[i] %= nc;
00840         }
00841 
00842         /* split offd part */
00843         hypre_CSRMatrixSplit( Aoffd, nr, nc, row_block_num, offd_col_block_num, csr_blocks );
00844 
00845         for( i = 0; i < num_blocks; i++ )
00846         {
00847             hypre_TFree( hypre_ParCSRMatrixOffd( blocks[i] ) );
00848             hypre_ParCSRMatrixOffd( blocks[i] ) = csr_blocks[i];
00849         }
00850 
00851         hypre_TFree( csr_blocks );
00852         hypre_TFree( col_block_num );
00853         hypre_TFree( row_block_num );
00854 
00855         /* update block col-maps */
00856         for( int bi = 0; bi < nr; bi++ )
00857         {
00858             for( int bj = 0; bj < nc; bj++ )
00859             {
00860                 hypre_ParCSRMatrix* block   = blocks[bi * nc + bj];
00861                 hypre_CSRMatrix* block_offd = hypre_ParCSRMatrixOffd( block );
00862                 HYPRE_Int block_offd_cols   = hypre_CSRMatrixNumCols( block_offd );
00863 
00864                 HYPRE_Int* block_col_map = hypre_TAlloc( HYPRE_Int, block_offd_cols );
00865                 for( i = j = 0; i < offd_cols; i++ )
00866                 {
00867                     HYPRE_Int bn = offd_col_block_num[i];
00868                     if( bn == bj )
00869                     {
00870                         block_col_map[j++] = offd_global_col[i];
00871                     }
00872                 }
00873                 hypre_assert( j == block_offd_cols );
00874 
00875                 hypre_ParCSRMatrixColMapOffd( block ) = block_col_map;
00876             }
00877         }
00878 
00879         hypre_TFree( offd_global_col );
00880         hypre_TFree( offd_col_block_num );
00881 
00882         /* finish the new matrices, make them own all the stuff */
00883         for( i = 0; i < num_blocks; i++ )
00884         {
00885             hypre_ParCSRMatrixSetNumNonzeros( blocks[i] );
00886             hypre_MatvecCommPkgCreate( blocks[i] );
00887 
00888             hypre_ParCSRMatrixOwnsData( blocks[i] ) = 1;
00889 
00890             /* only the first block will own the row/col_starts */
00891             hypre_ParCSRMatrixOwnsRowStarts( blocks[i] ) = !i;
00892             hypre_ParCSRMatrixOwnsColStarts( blocks[i] ) = !i;
00893         }
00894     }
00895 
00896     /* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
00897     void hypre_CSRMatrixBooleanMatvec( hypre_CSRMatrix* A,
00898                                        HYPRE_Bool alpha,
00899                                        HYPRE_Bool* x,
00900                                        HYPRE_Bool beta,
00901                                        HYPRE_Bool* y )
00902     {
00903         /* HYPRE_Complex    *A_data   = hypre_CSRMatrixData(A); */
00904         HYPRE_Int* A_i     = hypre_CSRMatrixI( A );
00905         HYPRE_Int* A_j     = hypre_CSRMatrixJ( A );
00906         HYPRE_Int num_rows = hypre_CSRMatrixNumRows( A );
00907 
00908         HYPRE_Int* A_rownnz  = hypre_CSRMatrixRownnz( A );
00909         HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz( A );
00910 
00911         HYPRE_Bool* x_data = x;
00912         HYPRE_Bool* y_data = y;
00913 
00914         HYPRE_Bool temp, tempx;
00915 
00916         HYPRE_Int i, jj;
00917 
00918         HYPRE_Int m;
00919 
00920         HYPRE_Real xpar = 0.7;
00921 
00922         /*-----------------------------------------------------------------------
00923          * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
00924          *-----------------------------------------------------------------------*/
00925 
00926         if( alpha == 0 )
00927         {
00928 #ifdef HYPRE_USING_OPENMP
00929 #pragma omp parallel for private( i ) HYPRE_SMP_SCHEDULE
00930 #endif
00931             for( i = 0; i < num_rows; i++ )
00932             {
00933                 y_data[i] = y_data[i] && beta;
00934             }
00935             return;
00936         }
00937 
00938         /*-----------------------------------------------------------------------
00939          * y = (beta/alpha)*y
00940          *-----------------------------------------------------------------------*/
00941 
00942         if( beta == 0 )
00943         {
00944 #ifdef HYPRE_USING_OPENMP
00945 #pragma omp parallel for private( i ) HYPRE_SMP_SCHEDULE
00946 #endif
00947             for( i = 0; i < num_rows; i++ )
00948             {
00949                 y_data[i] = 0;
00950             }
00951         }
00952         else
00953         {
00954             /* beta is true -> no change to y_data */
00955         }
00956 
00957         /*-----------------------------------------------------------------
00958          * y += A*x
00959          *-----------------------------------------------------------------*/
00960 
00961         /* use rownnz pointer to do the A*x multiplication  when num_rownnz is smaller than num_rows
00962          */
00963 
00964         if( num_rownnz < xpar * ( num_rows ) )
00965         {
00966 #ifdef HYPRE_USING_OPENMP
00967 #pragma omp parallel for private( i, jj, m, tempx ) HYPRE_SMP_SCHEDULE
00968 #endif
00969             for( i = 0; i < num_rownnz; i++ )
00970             {
00971                 m = A_rownnz[i];
00972 
00973                 tempx = 0;
00974                 for( jj = A_i[m]; jj < A_i[m + 1]; jj++ )
00975                 {
00976                     /* tempx = tempx || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
00977                     tempx = tempx || x_data[A_j[jj]];
00978                 }
00979                 y_data[m] = y_data[m] || tempx;
00980             }
00981         }
00982         else
00983         {
00984 #ifdef HYPRE_USING_OPENMP
00985 #pragma omp parallel for private( i, jj, temp ) HYPRE_SMP_SCHEDULE
00986 #endif
00987             for( i = 0; i < num_rows; i++ )
00988             {
00989                 temp = 0;
00990                 for( jj = A_i[i]; jj < A_i[i + 1]; jj++ )
00991                 {
00992                     /* temp = temp || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
00993                     temp = temp || x_data[A_j[jj]];
00994                 }
00995                 y_data[i] = y_data[i] || temp;
00996             }
00997         }
00998 
00999         /*-----------------------------------------------------------------
01000          * y = alpha*y
01001          *-----------------------------------------------------------------*/
01002         /* alpha is true */
01003     }
01004 
01005     /* Based on hypre_ParCSRCommHandleCreate in hypre's par_csr_communication.c. The
01006         input variable job controls the communication type: 1=Matvec, 2=MatvecT. */
01007     hypre_ParCSRCommHandle* hypre_ParCSRCommHandleCreate_bool( HYPRE_Int job,
01008                                                                hypre_ParCSRCommPkg* comm_pkg,
01009                                                                HYPRE_Bool* send_data,
01010                                                                HYPRE_Bool* recv_data )
01011     {
01012         HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
01013         HYPRE_Int num_recvs = hypre_ParCSRCommPkgNumRecvs( comm_pkg );
01014         MPI_Comm comm       = hypre_ParCSRCommPkgComm( comm_pkg );
01015 
01016         hypre_ParCSRCommHandle* comm_handle;
01017         HYPRE_Int num_requests;
01018         hypre_MPI_Request* requests;
01019 
01020         HYPRE_Int i, j;
01021         HYPRE_Int my_id, num_procs;
01022         HYPRE_Int ip, vec_start, vec_len;
01023 
01024         num_requests = num_sends + num_recvs;
01025         requests     = hypre_CTAlloc( hypre_MPI_Request, num_requests );
01026 
01027         hypre_MPI_Comm_size( comm, &num_procs );
01028         hypre_MPI_Comm_rank( comm, &my_id );
01029 
01030         j = 0;
01031         switch( job )
01032         {
01033             case 1: {
01034                 HYPRE_Bool* d_send_data = (HYPRE_Bool*)send_data;
01035                 HYPRE_Bool* d_recv_data = (HYPRE_Bool*)recv_data;
01036                 for( i = 0; i < num_recvs; i++ )
01037                 {
01038                     ip        = hypre_ParCSRCommPkgRecvProc( comm_pkg, i );
01039                     vec_start = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i );
01040                     vec_len   = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i + 1 ) - vec_start;
01041                     hypre_MPI_Irecv( &d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
01042                 }
01043                 for( i = 0; i < num_sends; i++ )
01044                 {
01045                     vec_start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
01046                     vec_len   = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ) - vec_start;
01047                     ip        = hypre_ParCSRCommPkgSendProc( comm_pkg, i );
01048                     hypre_MPI_Isend( &d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
01049                 }
01050                 break;
01051             }
01052             case 2: {
01053                 HYPRE_Bool* d_send_data = (HYPRE_Bool*)send_data;
01054                 HYPRE_Bool* d_recv_data = (HYPRE_Bool*)recv_data;
01055                 for( i = 0; i < num_sends; i++ )
01056                 {
01057                     vec_start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
01058                     vec_len   = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ) - vec_start;
01059                     ip        = hypre_ParCSRCommPkgSendProc( comm_pkg, i );
01060                     hypre_MPI_Irecv( &d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
01061                 }
01062                 for( i = 0; i < num_recvs; i++ )
01063                 {
01064                     ip        = hypre_ParCSRCommPkgRecvProc( comm_pkg, i );
01065                     vec_start = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i );
01066                     vec_len   = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i + 1 ) - vec_start;
01067                     hypre_MPI_Isend( &d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
01068                 }
01069                 break;
01070             }
01071         }
01072         /*--------------------------------------------------------------------
01073          * set up comm_handle and return
01074          *--------------------------------------------------------------------*/
01075 
01076         comm_handle = hypre_CTAlloc( hypre_ParCSRCommHandle, 1 );
01077 
01078         hypre_ParCSRCommHandleCommPkg( comm_handle )     = comm_pkg;
01079         hypre_ParCSRCommHandleSendData( comm_handle )    = send_data;
01080         hypre_ParCSRCommHandleRecvData( comm_handle )    = recv_data;
01081         hypre_ParCSRCommHandleNumRequests( comm_handle ) = num_requests;
01082         hypre_ParCSRCommHandleRequests( comm_handle )    = requests;
01083 
01084         return comm_handle;
01085     }
01086 
01087     /* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
01088     void hypre_ParCSRMatrixBooleanMatvec( hypre_ParCSRMatrix* A,
01089                                           HYPRE_Bool alpha,
01090                                           HYPRE_Bool* x,
01091                                           HYPRE_Bool beta,
01092                                           HYPRE_Bool* y )
01093     {
01094         hypre_ParCSRCommHandle* comm_handle;
01095         hypre_ParCSRCommPkg* comm_pkg = hypre_ParCSRMatrixCommPkg( A );
01096         hypre_CSRMatrix* diag         = hypre_ParCSRMatrixDiag( A );
01097         hypre_CSRMatrix* offd         = hypre_ParCSRMatrixOffd( A );
01098 
01099         HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols( offd );
01100         HYPRE_Int num_sends, i, j, index;
01101 
01102         HYPRE_Bool *x_tmp, *x_buf;
01103 
01104         x_tmp = hypre_CTAlloc( HYPRE_Bool, num_cols_offd );
01105 
01106         /*---------------------------------------------------------------------
01107          * If there exists no CommPkg for A, a CommPkg is generated using
01108          * equally load balanced partitionings
01109          *--------------------------------------------------------------------*/
01110         if( !comm_pkg )
01111         {
01112             hypre_MatvecCommPkgCreate( A );
01113             comm_pkg = hypre_ParCSRMatrixCommPkg( A );
01114         }
01115 
01116         num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
01117         x_buf     = hypre_CTAlloc( HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
01118 
01119         index = 0;
01120         for( i = 0; i < num_sends; i++ )
01121         {
01122             j = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
01123             for( ; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
01124             {
01125                 x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j )];
01126             }
01127         }
01128 
01129         comm_handle = hypre_ParCSRCommHandleCreate_bool( 1, comm_pkg, x_buf, x_tmp );
01130 
01131         hypre_CSRMatrixBooleanMatvec( diag, alpha, x, beta, y );
01132 
01133         hypre_ParCSRCommHandleDestroy( comm_handle );
01134 
01135         if( num_cols_offd )
01136         {
01137             hypre_CSRMatrixBooleanMatvec( offd, alpha, x_tmp, 1, y );
01138         }
01139 
01140         hypre_TFree( x_buf );
01141         hypre_TFree( x_tmp );
01142     }
01143 
01144     HYPRE_Int hypre_CSRMatrixSum( hypre_CSRMatrix* A, HYPRE_Complex beta, hypre_CSRMatrix* B )
01145     {
01146         HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
01147         HYPRE_Int* A_i        = hypre_CSRMatrixI( A );
01148         HYPRE_Int* A_j        = hypre_CSRMatrixJ( A );
01149         HYPRE_Int nrows_A     = hypre_CSRMatrixNumRows( A );
01150         HYPRE_Int ncols_A     = hypre_CSRMatrixNumCols( A );
01151         HYPRE_Complex* B_data = hypre_CSRMatrixData( B );
01152         HYPRE_Int* B_i        = hypre_CSRMatrixI( B );
01153         HYPRE_Int* B_j        = hypre_CSRMatrixJ( B );
01154         HYPRE_Int nrows_B     = hypre_CSRMatrixNumRows( B );
01155         HYPRE_Int ncols_B     = hypre_CSRMatrixNumCols( B );
01156 
01157         HYPRE_Int ia, j, pos;
01158         HYPRE_Int* marker;
01159 
01160         if( nrows_A != nrows_B || ncols_A != ncols_B )
01161         {
01162             return -1; /* error: incompatible matrix dimensions */
01163         }
01164 
01165         marker = hypre_CTAlloc( HYPRE_Int, ncols_A );
01166         for( ia = 0; ia < ncols_A; ia++ )
01167         {
01168             marker[ia] = -1;
01169         }
01170 
01171         for( ia = 0; ia < nrows_A; ia++ )
01172         {
01173             for( j = A_i[ia]; j < A_i[ia + 1]; j++ )
01174             {
01175                 marker[A_j[j]] = j;
01176             }
01177 
01178             for( j = B_i[ia]; j < B_i[ia + 1]; j++ )
01179             {
01180                 pos = marker[B_j[j]];
01181                 if( pos < A_i[ia] )
01182                 {
01183                     return -2; /* error: found an entry in B that is not present in A */
01184                 }
01185                 A_data[pos] += beta * B_data[j];
01186             }
01187         }
01188 
01189         hypre_TFree( marker );
01190         return 0;
01191     }
01192 
01193     hypre_ParCSRMatrix* hypre_ParCSRMatrixAdd( hypre_ParCSRMatrix* A, hypre_ParCSRMatrix* B )
01194     {
01195         MPI_Comm comm           = hypre_ParCSRMatrixComm( A );
01196         hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
01197         hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
01198         HYPRE_Int* A_cmap       = hypre_ParCSRMatrixColMapOffd( A );
01199         HYPRE_Int A_cmap_size   = hypre_CSRMatrixNumCols( A_offd );
01200         hypre_CSRMatrix* B_diag = hypre_ParCSRMatrixDiag( B );
01201         hypre_CSRMatrix* B_offd = hypre_ParCSRMatrixOffd( B );
01202         HYPRE_Int* B_cmap       = hypre_ParCSRMatrixColMapOffd( B );
01203         HYPRE_Int B_cmap_size   = hypre_CSRMatrixNumCols( B_offd );
01204         hypre_ParCSRMatrix* C;
01205         hypre_CSRMatrix* C_diag;
01206         hypre_CSRMatrix* C_offd;
01207         HYPRE_Int* C_cmap;
01208         HYPRE_Int im;
01209         HYPRE_Int cmap_differ;
01210 
01211         /* Check if A_cmap and B_cmap are the same. */
01212         cmap_differ = 0;
01213         if( A_cmap_size != B_cmap_size )
01214         {
01215             cmap_differ = 1; /* A and B have different cmap_size */
01216         }
01217         else
01218         {
01219             for( im = 0; im < A_cmap_size; im++ )
01220             {
01221                 if( A_cmap[im] != B_cmap[im] )
01222                 {
01223                     cmap_differ = 1; /* A and B have different cmap arrays */
01224                     break;
01225                 }
01226             }
01227         }
01228 
01229         if( cmap_differ == 0 )
01230         {
01231             /* A and B have the same column mapping for their off-diagonal blocks so
01232                 we can sum the diagonal and off-diagonal blocks separately and reduce
01233                 temporary memory usage. */
01234 
01235             /* Add diagonals, off-diagonals, copy cmap. */
01236             C_diag = hypre_CSRMatrixAdd( A_diag, B_diag );
01237             if( !C_diag )
01238             {
01239                 return NULL; /* error: A_diag and B_diag have different dimensions */
01240             }
01241             C_offd = hypre_CSRMatrixAdd( A_offd, B_offd );
01242             if( !C_offd )
01243             {
01244                 hypre_CSRMatrixDestroy( C_diag );
01245                 return NULL; /* error: A_offd and B_offd have different dimensions */
01246             }
01247             /* copy A_cmap -> C_cmap */
01248             C_cmap = hypre_TAlloc( HYPRE_Int, A_cmap_size );
01249             for( im = 0; im < A_cmap_size; im++ )
01250             {
01251                 C_cmap[im] = A_cmap[im];
01252             }
01253 
01254             C = hypre_ParCSRMatrixCreate( comm, hypre_ParCSRMatrixGlobalNumRows( A ),
01255                                           hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
01256                                           hypre_ParCSRMatrixColStarts( A ), hypre_CSRMatrixNumCols( C_offd ),
01257                                           hypre_CSRMatrixNumNonzeros( C_diag ), hypre_CSRMatrixNumNonzeros( C_offd ) );
01258 
01259             /* In C, destroy diag/offd (allocated by Create) and replace them with
01260             C_diag/C_offd. */
01261             hypre_CSRMatrixDestroy( hypre_ParCSRMatrixDiag( C ) );
01262             hypre_CSRMatrixDestroy( hypre_ParCSRMatrixOffd( C ) );
01263             hypre_ParCSRMatrixDiag( C ) = C_diag;
01264             hypre_ParCSRMatrixOffd( C ) = C_offd;
01265 
01266             hypre_ParCSRMatrixColMapOffd( C ) = C_cmap;
01267         }
01268         else
01269         {
01270             /* A and B have different column mappings for their off-diagonal blocks so
01271             we need to use the column maps to create full-width CSR matricies. */
01272 
01273             int ierr = 0;
01274             hypre_CSRMatrix* csr_A;
01275             hypre_CSRMatrix* csr_B;
01276             hypre_CSRMatrix* csr_C_temp;
01277 
01278             /* merge diag and off-diag portions of A */
01279             csr_A = hypre_MergeDiagAndOffd( A );
01280 
01281             /* merge diag and off-diag portions of B */
01282             csr_B = hypre_MergeDiagAndOffd( B );
01283 
01284             /* add A and B */
01285             csr_C_temp = hypre_CSRMatrixAdd( csr_A, csr_B );
01286 
01287             /* delete CSR versions of A and B */
01288             ierr += hypre_CSRMatrixDestroy( csr_A );
01289             ierr += hypre_CSRMatrixDestroy( csr_B );
01290 
01291             /* create a new empty ParCSR matrix to contain the sum */
01292             C = hypre_ParCSRMatrixCreate( hypre_ParCSRMatrixComm( A ), hypre_ParCSRMatrixGlobalNumRows( A ),
01293                                           hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
01294                                           hypre_ParCSRMatrixColStarts( A ), 0, 0, 0 );
01295 
01296             /* split C into diag and off-diag portions */
01297             /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the
01298                 number of columns in csr_C_temp which is the global number of columns
01299                 in A and B. This does not scale well. */
01300             ierr += GenerateDiagAndOffd( csr_C_temp, C, hypre_ParCSRMatrixFirstColDiag( A ),
01301                                          hypre_ParCSRMatrixLastColDiag( A ) );
01302 
01303             /* delete CSR version of C */
01304             ierr += hypre_CSRMatrixDestroy( csr_C_temp );
01305         }
01306 
01307         /* hypre_ParCSRMatrixSetNumNonzeros(A); */
01308 
01309         /* Make sure that the first entry in each row is the diagonal one. */
01310         hypre_CSRMatrixReorder( hypre_ParCSRMatrixDiag( C ) );
01311 
01312         /* C owns diag, offd, and cmap. */
01313         hypre_ParCSRMatrixSetDataOwner( C, 1 );
01314         /* C does not own row and column starts. */
01315         hypre_ParCSRMatrixSetRowStartsOwner( C, 0 );
01316         hypre_ParCSRMatrixSetColStartsOwner( C, 0 );
01317 
01318         return C;
01319     }
01320 
01321     HYPRE_Int hypre_ParCSRMatrixSum( hypre_ParCSRMatrix* A, HYPRE_Complex beta, hypre_ParCSRMatrix* B )
01322     {
01323         hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
01324         hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
01325         hypre_CSRMatrix* B_diag = hypre_ParCSRMatrixDiag( B );
01326         hypre_CSRMatrix* B_offd = hypre_ParCSRMatrixOffd( B );
01327         HYPRE_Int error;
01328 
01329         error = hypre_CSRMatrixSum( A_diag, beta, B_diag );
01330         error = error ? error : hypre_CSRMatrixSum( A_offd, beta, B_offd );
01331 
01332         return error;
01333     }
01334 
01335     HYPRE_Int hypre_CSRMatrixSetConstantValues( hypre_CSRMatrix* A, HYPRE_Complex value )
01336     {
01337         HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
01338         HYPRE_Int A_nnz       = hypre_CSRMatrixNumNonzeros( A );
01339         HYPRE_Int ia;
01340 
01341         for( ia = 0; ia < A_nnz; ia++ )
01342         {
01343             A_data[ia] = value;
01344         }
01345 
01346         return 0;
01347     }
01348 
01349     HYPRE_Int hypre_ParCSRMatrixSetConstantValues( hypre_ParCSRMatrix* A, HYPRE_Complex value )
01350     {
01351         hypre_CSRMatrixSetConstantValues( hypre_ParCSRMatrixDiag( A ), value );
01352         hypre_CSRMatrixSetConstantValues( hypre_ParCSRMatrixOffd( A ), value );
01353 
01354         return 0;
01355     }
01356 
01357 }  // namespace internal
01358 
01359 }  // namespace moab
01360 
01361 #endif  // MOAB_HAVE_MPI
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines