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