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