 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
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
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
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.