 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11
12 #include "hypre_parcsr.hpp"
13
14 #ifdef MOAB_HAVE_MPI
15
16 #include <limits>
17
18 namespace moab
19 {
20
21 namespace internal
22 {
23
24 /*--------------------------------------------------------------------------
25 * A*X = B style elimination
26 *--------------------------------------------------------------------------*/
27
28 /*
29 Function: hypre_CSRMatrixEliminateAXB
30
31 Eliminate the rows and columns of A corresponding to the
32 given sorted (!) list of rows. Put I on the eliminated diagonal,
33 subtract columns times X from B.
34 */
35 void hypre_CSRMatrixEliminateAXB( hypre_CSRMatrix* A,
36 HYPRE_Int nrows_to_eliminate,
37 HYPRE_Int* rows_to_eliminate,
38 hypre_Vector* X,
39 hypre_Vector* B )
40 {
41 HYPRE_Int i, j;
42 HYPRE_Int irow, jcol, ibeg, iend, pos;
43 HYPRE_Real a;
44
45 HYPRE_Int* Ai = hypre_CSRMatrixI( A );
46 HYPRE_Int* Aj = hypre_CSRMatrixJ( A );
47 HYPRE_Real* Adata = hypre_CSRMatrixData( A );
48 HYPRE_Int nrows = hypre_CSRMatrixNumRows( A );
49
50 HYPRE_Real* Xdata = hypre_VectorData( X );
51 HYPRE_Real* Bdata = hypre_VectorData( B );
52
53 /* eliminate the columns */
54 for( i = 0; i < nrows; i++ )
55 {
56 ibeg = Ai[i];
57 iend = Ai[i + 1];
58 for( j = ibeg; j < iend; j++ )
59 {
60 jcol = Aj[j];
61 pos = hypre_BinarySearch( rows_to_eliminate, jcol, nrows_to_eliminate );
62 if( pos != -1 )
63 {
64 a = Adata[j];
65 Adata[j] = 0.0;
66 Bdata[i] -= a * Xdata[jcol];
67 }
68 }
69 }
70
71 /* remove the rows and set the diagonal equal to 1 */
72 for( i = 0; i < nrows_to_eliminate; i++ )
73 {
74 irow = rows_to_eliminate[i];
75 ibeg = Ai[irow];
76 iend = Ai[irow + 1];
77 for( j = ibeg; j < iend; j++ )
78 {
79 if( Aj[j] == irow )
80 {
81 Adata[j] = 1.0;
82 }
83 else
84 {
85 Adata[j] = 0.0;
86 }
87 }
88 }
89 }
90
91 /*
92 Function: hypre_CSRMatrixEliminateOffdColsAXB
93
94 Eliminate the given sorted (!) list of columns of A, subtract them from B.
95 */
96 void hypre_CSRMatrixEliminateOffdColsAXB( hypre_CSRMatrix* A,
97 HYPRE_Int ncols_to_eliminate,
98 HYPRE_Int* eliminate_cols,
99 HYPRE_Real* eliminate_coefs,
100 hypre_Vector* B )
101 {
102 HYPRE_Int i, j;
103 HYPRE_Int ibeg, iend, pos;
104 HYPRE_Real a;
105
106 HYPRE_Int* Ai = hypre_CSRMatrixI( A );
107 HYPRE_Int* Aj = hypre_CSRMatrixJ( A );
108 HYPRE_Real* Adata = hypre_CSRMatrixData( A );
109 HYPRE_Int nrows = hypre_CSRMatrixNumRows( A );
110
111 HYPRE_Real* Bdata = hypre_VectorData( B );
112
113 for( i = 0; i < nrows; i++ )
114 {
115 ibeg = Ai[i];
116 iend = Ai[i + 1];
117 for( j = ibeg; j < iend; j++ )
118 {
119 pos = hypre_BinarySearch( eliminate_cols, Aj[j], ncols_to_eliminate );
120 if( pos != -1 )
121 {
122 a = Adata[j];
123 Adata[j] = 0.0;
124 Bdata[i] -= a * eliminate_coefs[pos];
125 }
126 }
127 }
128 }
129
130 /*
131 Function: hypre_CSRMatrixEliminateOffdRowsAXB
132
133 Eliminate (zero) the given list of rows of A.
134 */
135 void hypre_CSRMatrixEliminateOffdRowsAXB( hypre_CSRMatrix* A,
136 HYPRE_Int nrows_to_eliminate,
137 HYPRE_Int* rows_to_eliminate )
138 {
139 HYPRE_Int* Ai = hypre_CSRMatrixI( A );
140 HYPRE_Real* Adata = hypre_CSRMatrixData( A );
141
142 HYPRE_Int i, j;
143 HYPRE_Int irow, ibeg, iend;
144
145 for( i = 0; i < nrows_to_eliminate; i++ )
146 {
147 irow = rows_to_eliminate[i];
148 ibeg = Ai[irow];
149 iend = Ai[irow + 1];
150 for( j = ibeg; j < iend; j++ )
151 {
152 Adata[j] = 0.0;
153 }
154 }
155 }
156
157 /*
158 Function: hypre_ParCSRMatrixEliminateAXB
159
160 This function eliminates the global rows and columns of a matrix
161 A corresponding to given lists of sorted (!) local row numbers,
162 so that the solution to the system A*X = B is X_b for the given rows.
163
164 The elimination is done as follows:
165
166 (input) (output)
167
168 / A_ii | A_ib \ / A_ii | 0 \
169 A = | -----+----- | ---> | -----+----- |
170 \ A_bi | A_bb / \ 0 | I /
171
172 / X_i \ / X_i \
173 X = | --- | ---> | --- | (no change)
174 \ X_b / \ X_b /
175
176 / B_i \ / B_i - A_ib * X_b \
177 B = | --- | ---> | ---------------- |
178 \ B_b / \ X_b /
179
180 */
181 void hypre_ParCSRMatrixEliminateAXB( hypre_ParCSRMatrix* A,
182 HYPRE_Int num_rowscols_to_elim,
183 HYPRE_Int* rowscols_to_elim,
184 hypre_ParVector* X,
185 hypre_ParVector* B )
186 {
187 hypre_CSRMatrix* diag = hypre_ParCSRMatrixDiag( A );
188 hypre_CSRMatrix* offd = hypre_ParCSRMatrixOffd( A );
189 HYPRE_Int diag_nrows = hypre_CSRMatrixNumRows( diag );
190 HYPRE_Int offd_ncols = hypre_CSRMatrixNumCols( offd );
191
192 hypre_Vector* Xlocal = hypre_ParVectorLocalVector( X );
193 hypre_Vector* Blocal = hypre_ParVectorLocalVector( B );
194
195 HYPRE_Real* Bdata = hypre_VectorData( Blocal );
196 HYPRE_Real* Xdata = hypre_VectorData( Xlocal );
197
198 HYPRE_Int num_offd_cols_to_elim;
199 HYPRE_Int* offd_cols_to_elim;
200 HYPRE_Real* eliminate_coefs;
201
202 /* figure out which offd cols should be eliminated and with what coef */
203 hypre_ParCSRCommHandle* comm_handle;
204 hypre_ParCSRCommPkg* comm_pkg;
205 HYPRE_Int num_sends;
206 HYPRE_Int index, start;
207 HYPRE_Int i, j, k, irow;
208
209 HYPRE_Real* eliminate_row = hypre_CTAlloc( HYPRE_Real, diag_nrows );
210 HYPRE_Real* eliminate_col = hypre_CTAlloc( HYPRE_Real, offd_ncols );
211 HYPRE_Real *buf_data, coef;
212
213 /* make sure A has a communication package */
214 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
215 if( !comm_pkg )
216 {
217 hypre_MatvecCommPkgCreate( A );
218 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
219 }
220
221 /* HACK: rows that shouldn't be eliminated are marked with quiet NaN;
222 those that should are set to the boundary value from X; this is to
223 avoid sending complex type (int+double) or communicating twice. */
224 for( i = 0; i < diag_nrows; i++ )
225 {
226 eliminate_row[i] = std::numeric_limits< HYPRE_Real >::quiet_NaN();
227 }
228 for( i = 0; i < num_rowscols_to_elim; i++ )
229 {
230 irow = rowscols_to_elim[i];
231 eliminate_row[irow] = Xdata[irow];
232 }
233
234 /* use a Matvec communication pattern to find (in eliminate_col)
235 which of the local offd columns are to be eliminated */
236 num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
237 buf_data = hypre_CTAlloc( HYPRE_Real, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
238 index = 0;
239 for( i = 0; i < num_sends; i++ )
240 {
241 start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
242 for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
243 {
244 k = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
245 buf_data[index++] = eliminate_row[k];
246 }
247 }
248 comm_handle = hypre_ParCSRCommHandleCreate( 1, comm_pkg, buf_data, eliminate_col );
249
250 /* do sequential part of the elimination while stuff is getting sent */
251 hypre_CSRMatrixEliminateAXB( diag, num_rowscols_to_elim, rowscols_to_elim, Xlocal, Blocal );
252
253 /* finish the communication */
254 hypre_ParCSRCommHandleDestroy( comm_handle );
255
256 /* received eliminate_col[], count offd columns to eliminate */
257 num_offd_cols_to_elim = 0;
258 for( i = 0; i < offd_ncols; i++ )
259 {
260 coef = eliminate_col[i];
261 if( coef == coef ) // test for NaN
262 {
263 num_offd_cols_to_elim++;
264 }
265 }
266
267 offd_cols_to_elim = hypre_CTAlloc( HYPRE_Int, num_offd_cols_to_elim );
268 eliminate_coefs = hypre_CTAlloc( HYPRE_Real, num_offd_cols_to_elim );
269
270 /* get a list of offd column indices and coefs */
271 num_offd_cols_to_elim = 0;
272 for( i = 0; i < offd_ncols; i++ )
273 {
274 coef = eliminate_col[i];
275 if( coef == coef ) // test for NaN
276 {
277 offd_cols_to_elim[num_offd_cols_to_elim] = i;
278 eliminate_coefs[num_offd_cols_to_elim] = coef;
279 num_offd_cols_to_elim++;
280 }
281 }
282
283 hypre_TFree( buf_data );
284 hypre_TFree( eliminate_row );
285 hypre_TFree( eliminate_col );
286
287 /* eliminate the off-diagonal part */
288 hypre_CSRMatrixEliminateOffdColsAXB( offd, num_offd_cols_to_elim, offd_cols_to_elim, eliminate_coefs, Blocal );
289
290 hypre_CSRMatrixEliminateOffdRowsAXB( offd, num_rowscols_to_elim, rowscols_to_elim );
291
292 /* set boundary values in the rhs */
293 for( int m = 0; m < num_rowscols_to_elim; m++ )
294 {
295 irow = rowscols_to_elim[m];
296 Bdata[irow] = Xdata[irow];
297 }
298
299 hypre_TFree( offd_cols_to_elim );
300 hypre_TFree( eliminate_coefs );
301 }
302
303 /*--------------------------------------------------------------------------
304 * (A + Ae) style elimination
305 *--------------------------------------------------------------------------*/
306
307 /*
308 Function: hypre_CSRMatrixElimCreate
309
310 Prepare the Ae matrix: count nnz, initialize I, allocate J and data.
311 */
312 void hypre_CSRMatrixElimCreate( hypre_CSRMatrix* A,
313 hypre_CSRMatrix* Ae,
314 HYPRE_Int nrows,
315 HYPRE_Int* rows,
316 HYPRE_Int ncols,
317 HYPRE_Int* cols,
318 HYPRE_Int* col_mark )
319 {
320 HYPRE_Int i, j, col;
321 HYPRE_Int A_beg, A_end;
322
323 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
324 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
325 HYPRE_Int A_rows = hypre_CSRMatrixNumRows( A );
326
327 hypre_CSRMatrixI( Ae ) = hypre_TAlloc( HYPRE_Int, A_rows + 1 );
328
329 HYPRE_Int* Ae_i = hypre_CSRMatrixI( Ae );
330 HYPRE_Int nnz = 0;
331
332 for( i = 0; i < A_rows; i++ )
333 {
334 Ae_i[i] = nnz;
335
336 A_beg = A_i[i];
337 A_end = A_i[i + 1];
338
339 if( hypre_BinarySearch( rows, i, nrows ) >= 0 )
340 {
341 /* full row */
342 nnz += A_end - A_beg;
343
344 if( col_mark )
345 {
346 for( j = A_beg; j < A_end; j++ )
347 {
348 col_mark[A_j[j]] = 1;
349 }
350 }
351 }
352 else
353 {
354 /* count columns */
355 for( j = A_beg; j < A_end; j++ )
356 {
357 col = A_j[j];
358 if( hypre_BinarySearch( cols, col, ncols ) >= 0 )
359 {
360 nnz++;
361 if( col_mark )
362 {
363 col_mark[col] = 1;
364 }
365 }
366 }
367 }
368 }
369 Ae_i[A_rows] = nnz;
370
371 hypre_CSRMatrixJ( Ae ) = hypre_TAlloc( HYPRE_Int, nnz );
372 hypre_CSRMatrixData( Ae ) = hypre_TAlloc( HYPRE_Real, nnz );
373 hypre_CSRMatrixNumNonzeros( Ae ) = nnz;
374 }
375
376 /*
377 Function: hypre_CSRMatrixEliminateRowsCols
378
379 Eliminate rows and columns of A, store eliminated values in Ae.
380 If 'diag' is nonzero, the eliminated diagonal of A is set to identity.
381 If 'col_remap' is not NULL it specifies renumbering of columns of Ae.
382 */
383 void hypre_CSRMatrixEliminateRowsCols( hypre_CSRMatrix* A,
384 hypre_CSRMatrix* Ae,
385 HYPRE_Int nrows,
386 HYPRE_Int* rows,
387 HYPRE_Int ncols,
388 HYPRE_Int* cols,
389 int diag,
390 HYPRE_Int* col_remap )
391 {
392 HYPRE_Int i, j, k, col;
393 HYPRE_Int A_beg, Ae_beg, A_end;
394 HYPRE_Real a;
395
396 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
397 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
398 HYPRE_Real* A_data = hypre_CSRMatrixData( A );
399 HYPRE_Int A_rows = hypre_CSRMatrixNumRows( A );
400
401 HYPRE_Int* Ae_i = hypre_CSRMatrixI( Ae );
402 HYPRE_Int* Ae_j = hypre_CSRMatrixJ( Ae );
403 HYPRE_Real* Ae_data = hypre_CSRMatrixData( Ae );
404
405 for( i = 0; i < A_rows; i++ )
406 {
407 A_beg = A_i[i];
408 A_end = A_i[i + 1];
409 Ae_beg = Ae_i[i];
410
411 if( hypre_BinarySearch( rows, i, nrows ) >= 0 )
412 {
413 /* eliminate row */
414 for( j = A_beg, k = Ae_beg; j < A_end; j++, k++ )
415 {
416 col = A_j[j];
417 Ae_j[k] = col_remap ? col_remap[col] : col;
418 a = ( diag && col == i ) ? 1.0 : 0.0;
419 Ae_data[k] = A_data[j] - a;
420 A_data[j] = a;
421 }
422 }
423 else
424 {
425 /* eliminate columns */
426 for( j = A_beg, k = Ae_beg; j < A_end; j++ )
427 {
428 col = A_j[j];
429 if( hypre_BinarySearch( cols, col, ncols ) >= 0 )
430 {
431 Ae_j[k] = col_remap ? col_remap[col] : col;
432 Ae_data[k] = A_data[j];
433 A_data[j] = 0.0;
434 k++;
435 }
436 }
437 }
438 }
439 }
440
441 /*
442 Function: hypre_ParCSRMatrixEliminateAAe
443
444 (input) (output)
445
446 / A_ii | A_ib \ / A_ii | 0 \
447 A = | -----+----- | ---> | -----+----- |
448 \ A_bi | A_bb / \ 0 | I /
449
450
451 / 0 | A_ib
452 \ Ae = | -----+--------- | \ A_bi | A_bb - I /
453
454 */
455 void hypre_ParCSRMatrixEliminateAAe( hypre_ParCSRMatrix* A,
456 hypre_ParCSRMatrix** Ae,
457 HYPRE_Int num_rowscols_to_elim,
458 HYPRE_Int* rowscols_to_elim )
459 {
460 HYPRE_Int i, j, k;
461
462 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
463 hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
464 HYPRE_Int A_diag_nrows = hypre_CSRMatrixNumRows( A_diag );
465 HYPRE_Int A_offd_ncols = hypre_CSRMatrixNumCols( A_offd );
466
467 *Ae = hypre_ParCSRMatrixCreate( hypre_ParCSRMatrixComm( A ), hypre_ParCSRMatrixGlobalNumRows( A ),
468 hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
469 hypre_ParCSRMatrixColStarts( A ), 0, 0, 0 );
470
471 hypre_ParCSRMatrixSetRowStartsOwner( *Ae, 0 );
472 hypre_ParCSRMatrixSetColStartsOwner( *Ae, 0 );
473
474 hypre_CSRMatrix* Ae_diag = hypre_ParCSRMatrixDiag( *Ae );
475 hypre_CSRMatrix* Ae_offd = hypre_ParCSRMatrixOffd( *Ae );
476 HYPRE_Int Ae_offd_ncols;
477
478 HYPRE_Int num_offd_cols_to_elim;
479 HYPRE_Int* offd_cols_to_elim;
480
481 HYPRE_Int* A_col_map_offd = hypre_ParCSRMatrixColMapOffd( A );
482 HYPRE_Int* Ae_col_map_offd;
483
484 HYPRE_Int* col_mark;
485 HYPRE_Int* col_remap;
486
487 /* figure out which offd cols should be eliminated */
488 {
489 hypre_ParCSRCommHandle* comm_handle;
490 hypre_ParCSRCommPkg* comm_pkg;
491 HYPRE_Int num_sends, *int_buf_data;
492 HYPRE_Int index, start;
493
494 HYPRE_Int* eliminate_row = hypre_CTAlloc( HYPRE_Int, A_diag_nrows );
495 HYPRE_Int* eliminate_col = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
496
497 /* make sure A has a communication package */
498 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
499 if( !comm_pkg )
500 {
501 hypre_MatvecCommPkgCreate( A );
502 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
503 }
504
505 /* which of the local rows are to be eliminated */
506 for( i = 0; i < A_diag_nrows; i++ )
507 {
508 eliminate_row[i] = 0;
509 }
510 for( i = 0; i < num_rowscols_to_elim; i++ )
511 {
512 eliminate_row[rowscols_to_elim[i]] = 1;
513 }
514
515 /* use a Matvec communication pattern to find (in eliminate_col)
516 which of the local offd columns are to be eliminated */
517 num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
518 int_buf_data = hypre_CTAlloc( HYPRE_Int, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
519 index = 0;
520 for( i = 0; i < num_sends; i++ )
521 {
522 start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
523 for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
524 {
525 k = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
526 int_buf_data[index++] = eliminate_row[k];
527 }
528 }
529 comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data, eliminate_col );
530
531 /* eliminate diagonal part, overlapping it with communication */
532 hypre_CSRMatrixElimCreate( A_diag, Ae_diag, num_rowscols_to_elim, rowscols_to_elim, num_rowscols_to_elim,
533 rowscols_to_elim, NULL );
534
535 hypre_CSRMatrixEliminateRowsCols( A_diag, Ae_diag, num_rowscols_to_elim, rowscols_to_elim,
536 num_rowscols_to_elim, rowscols_to_elim, 1, NULL );
537 hypre_CSRMatrixReorder( Ae_diag );
538
539 /* finish the communication */
540 hypre_ParCSRCommHandleDestroy( comm_handle );
541
542 /* received eliminate_col[], count offd columns to eliminate */
543 num_offd_cols_to_elim = 0;
544 for( i = 0; i < A_offd_ncols; i++ )
545 {
546 if( eliminate_col[i] )
547 {
548 num_offd_cols_to_elim++;
549 }
550 }
551
552 offd_cols_to_elim = hypre_CTAlloc( HYPRE_Int, num_offd_cols_to_elim );
553
554 /* get a list of offd column indices and coefs */
555 num_offd_cols_to_elim = 0;
556 for( i = 0; i < A_offd_ncols; i++ )
557 {
558 if( eliminate_col[i] )
559 {
560 offd_cols_to_elim[num_offd_cols_to_elim++] = i;
561 }
562 }
563
564 hypre_TFree( int_buf_data );
565 hypre_TFree( eliminate_row );
566 hypre_TFree( eliminate_col );
567 }
568
569 /* eliminate the off-diagonal part */
570 col_mark = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
571 col_remap = hypre_CTAlloc( HYPRE_Int, A_offd_ncols );
572
573 hypre_CSRMatrixElimCreate( A_offd, Ae_offd, num_rowscols_to_elim, rowscols_to_elim, num_offd_cols_to_elim,
574 offd_cols_to_elim, col_mark );
575
576 for( i = k = 0; i < A_offd_ncols; i++ )
577 {
578 if( col_mark[i] )
579 {
580 col_remap[i] = k++;
581 }
582 }
583
584 hypre_CSRMatrixEliminateRowsCols( A_offd, Ae_offd, num_rowscols_to_elim, rowscols_to_elim,
585 num_offd_cols_to_elim, offd_cols_to_elim, 0, col_remap );
586
587 /* create col_map_offd for Ae */
588 Ae_offd_ncols = 0;
589 for( i = 0; i < A_offd_ncols; i++ )
590 {
591 if( col_mark[i] )
592 {
593 Ae_offd_ncols++;
594 }
595 }
596
597 Ae_col_map_offd = hypre_CTAlloc( HYPRE_Int, Ae_offd_ncols );
598
599 Ae_offd_ncols = 0;
600 for( i = 0; i < A_offd_ncols; i++ )
601 {
602 if( col_mark[i] )
603 {
604 Ae_col_map_offd[Ae_offd_ncols++] = A_col_map_offd[i];
605 }
606 }
607
608 hypre_ParCSRMatrixColMapOffd( *Ae ) = Ae_col_map_offd;
609 hypre_CSRMatrixNumCols( Ae_offd ) = Ae_offd_ncols;
610
611 hypre_TFree( col_remap );
612 hypre_TFree( col_mark );
613 hypre_TFree( offd_cols_to_elim );
614
615 hypre_ParCSRMatrixSetNumNonzeros( *Ae );
616 hypre_MatvecCommPkgCreate( *Ae );
617 }
618
619 /*--------------------------------------------------------------------------
620 * Split
621 *--------------------------------------------------------------------------*/
622
623 void hypre_CSRMatrixSplit( hypre_CSRMatrix* A,
624 HYPRE_Int nr,
625 HYPRE_Int nc,
626 HYPRE_Int* row_block_num,
627 HYPRE_Int* col_block_num,
628 hypre_CSRMatrix** blocks )
629 {
630 HYPRE_Int i, j, k, bi, bj;
631
632 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
633 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
634 HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
635
636 HYPRE_Int A_rows = hypre_CSRMatrixNumRows( A );
637 HYPRE_Int A_cols = hypre_CSRMatrixNumCols( A );
638
639 HYPRE_Int* num_rows = hypre_CTAlloc( HYPRE_Int, nr );
640 HYPRE_Int* num_cols = hypre_CTAlloc( HYPRE_Int, nc );
641
642 HYPRE_Int* block_row = hypre_TAlloc( HYPRE_Int, A_rows );
643 HYPRE_Int* block_col = hypre_TAlloc( HYPRE_Int, A_cols );
644
645 for( i = 0; i < A_rows; i++ )
646 {
647 block_row[i] = num_rows[row_block_num[i]]++;
648 }
649 for( j = 0; j < A_cols; j++ )
650 {
651 block_col[j] = num_cols[col_block_num[j]]++;
652 }
653
654 /* allocate the blocks */
655 for( i = 0; i < nr; i++ )
656 {
657 for( j = 0; j < nc; j++ )
658 {
659 hypre_CSRMatrix* B = hypre_CSRMatrixCreate( num_rows[i], num_cols[j], 0 );
660 hypre_CSRMatrixI( B ) = hypre_CTAlloc( HYPRE_Int, num_rows[i] + 1 );
661 blocks[i * nc + j] = B;
662 }
663 }
664
665 /* count block row nnz */
666 for( i = 0; i < A_rows; i++ )
667 {
668 bi = row_block_num[i];
669 for( j = A_i[i]; j < A_i[i + 1]; j++ )
670 {
671 bj = col_block_num[A_j[j]];
672 hypre_CSRMatrix* B = blocks[bi * nc + bj];
673 hypre_CSRMatrixI( B )[block_row[i] + 1]++;
674 }
675 }
676
677 /* count block nnz */
678 for( k = 0; k < nr * nc; k++ )
679 {
680 hypre_CSRMatrix* B = blocks[k];
681 HYPRE_Int* B_i = hypre_CSRMatrixI( B );
682
683 HYPRE_Int nnz = 0, rs;
684 for( int m = 1; m <= hypre_CSRMatrixNumRows( B ); m++ )
685 {
686 rs = B_i[m], B_i[m] = nnz, nnz += rs;
687 }
688
689 hypre_CSRMatrixJ( B ) = hypre_TAlloc( HYPRE_Int, nnz );
690 hypre_CSRMatrixData( B ) = hypre_TAlloc( HYPRE_Complex, nnz );
691 hypre_CSRMatrixNumNonzeros( B ) = nnz;
692 }
693
694 /* populate blocks */
695 for( i = 0; i < A_rows; i++ )
696 {
697 bi = row_block_num[i];
698 for( j = A_i[i]; j < A_i[i + 1]; j++ )
699 {
700 k = A_j[j];
701 bj = col_block_num[k];
702 hypre_CSRMatrix* B = blocks[bi * nc + bj];
703 HYPRE_Int* bii = hypre_CSRMatrixI( B ) + block_row[i] + 1;
704 hypre_CSRMatrixJ( B )[*bii] = block_col[k];
705 hypre_CSRMatrixData( B )[*bii] = A_data[j];
706 ( *bii )++;
707 }
708 }
709
710 hypre_TFree( block_col );
711 hypre_TFree( block_row );
712
713 hypre_TFree( num_cols );
714 hypre_TFree( num_rows );
715 }
716
717 void hypre_ParCSRMatrixSplit( hypre_ParCSRMatrix* A,
718 HYPRE_Int nr,
719 HYPRE_Int nc,
720 hypre_ParCSRMatrix** blocks,
721 int interleaved_rows,
722 int interleaved_cols )
723 {
724 HYPRE_Int i, j, k;
725
726 MPI_Comm comm = hypre_ParCSRMatrixComm( A );
727
728 hypre_CSRMatrix* Adiag = hypre_ParCSRMatrixDiag( A );
729 hypre_CSRMatrix* Aoffd = hypre_ParCSRMatrixOffd( A );
730
731 HYPRE_Int global_rows = hypre_ParCSRMatrixGlobalNumRows( A );
732 HYPRE_Int global_cols = hypre_ParCSRMatrixGlobalNumCols( A );
733
734 HYPRE_Int local_rows = hypre_CSRMatrixNumRows( Adiag );
735 HYPRE_Int local_cols = hypre_CSRMatrixNumCols( Adiag );
736 HYPRE_Int offd_cols = hypre_CSRMatrixNumCols( Aoffd );
737
738 hypre_assert( local_rows % nr == 0 && local_cols % nc == 0 );
739 hypre_assert( global_rows % nr == 0 && global_cols % nc == 0 );
740
741 HYPRE_Int block_rows = local_rows / nr;
742 HYPRE_Int block_cols = local_cols / nc;
743 HYPRE_Int num_blocks = nr * nc;
744
745 /* mark local rows and columns with block number */
746 HYPRE_Int* row_block_num = hypre_TAlloc( HYPRE_Int, local_rows );
747 HYPRE_Int* col_block_num = hypre_TAlloc( HYPRE_Int, local_cols );
748
749 for( i = 0; i < local_rows; i++ )
750 {
751 row_block_num[i] = interleaved_rows ? ( i % nr ) : ( i / block_rows );
752 }
753 for( i = 0; i < local_cols; i++ )
754 {
755 col_block_num[i] = interleaved_cols ? ( i % nc ) : ( i / block_cols );
756 }
757
758 /* determine the block numbers for offd columns */
759 HYPRE_Int* offd_col_block_num = hypre_TAlloc( HYPRE_Int, offd_cols );
760 hypre_ParCSRCommHandle* comm_handle;
761 HYPRE_Int* int_buf_data;
762 {
763 /* make sure A has a communication package */
764 hypre_ParCSRCommPkg* comm_pkg = hypre_ParCSRMatrixCommPkg( A );
765 if( !comm_pkg )
766 {
767 hypre_MatvecCommPkgCreate( A );
768 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
769 }
770
771 /* calculate the final global column numbers for each block */
772 HYPRE_Int* count = hypre_CTAlloc( HYPRE_Int, nc );
773 HYPRE_Int* block_global_col = hypre_TAlloc( HYPRE_Int, local_cols );
774 HYPRE_Int first_col = hypre_ParCSRMatrixFirstColDiag( A ) / nc;
775 for( i = 0; i < local_cols; i++ )
776 {
777 block_global_col[i] = first_col + count[col_block_num[i]]++;
778 }
779 hypre_TFree( count );
780
781 /* use a Matvec communication pattern to determine offd_col_block_num */
782 HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
783 int_buf_data = hypre_CTAlloc( HYPRE_Int, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
784 HYPRE_Int start, index = 0;
785 for( i = 0; i < num_sends; i++ )
786 {
787 start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
788 for( j = start; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
789 {
790 k = hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j );
791 int_buf_data[index++] = col_block_num[k] + nc * block_global_col[k];
792 }
793 }
794 hypre_TFree( block_global_col );
795
796 comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data, offd_col_block_num );
797 }
798
799 /* create the block matrices */
800 HYPRE_Int num_procs = 1;
801 if( !hypre_ParCSRMatrixAssumedPartition( A ) )
802 {
803 hypre_MPI_Comm_size( comm, &num_procs );
804 }
805
806 HYPRE_Int* row_starts = hypre_TAlloc( HYPRE_Int, num_procs + 1 );
807 HYPRE_Int* col_starts = hypre_TAlloc( HYPRE_Int, num_procs + 1 );
808 for( i = 0; i <= num_procs; i++ )
809 {
810 row_starts[i] = hypre_ParCSRMatrixRowStarts( A )[i] / nr;
811 col_starts[i] = hypre_ParCSRMatrixColStarts( A )[i] / nc;
812 }
813
814 for( i = 0; i < num_blocks; i++ )
815 {
816 blocks[i] =
817 hypre_ParCSRMatrixCreate( comm, global_rows / nr, global_cols / nc, row_starts, col_starts, 0, 0, 0 );
818 }
819
820 /* split diag part */
821 hypre_CSRMatrix** csr_blocks = hypre_TAlloc( hypre_CSRMatrix*, nr * nc );
822 hypre_CSRMatrixSplit( Adiag, nr, nc, row_block_num, col_block_num, csr_blocks );
823
824 for( i = 0; i < num_blocks; i++ )
825 {
826 hypre_TFree( hypre_ParCSRMatrixDiag( blocks[i] ) );
827 hypre_ParCSRMatrixDiag( blocks[i] ) = csr_blocks[i];
828 }
829
830 /* finish communication, receive offd_col_block_num */
831 hypre_ParCSRCommHandleDestroy( comm_handle );
832 hypre_TFree( int_buf_data );
833
834 /* decode global offd column numbers */
835 HYPRE_Int* offd_global_col = hypre_TAlloc( HYPRE_Int, offd_cols );
836 for( i = 0; i < offd_cols; i++ )
837 {
838 offd_global_col[i] = offd_col_block_num[i] / nc;
839 offd_col_block_num[i] %= nc;
840 }
841
842 /* split offd part */
843 hypre_CSRMatrixSplit( Aoffd, nr, nc, row_block_num, offd_col_block_num, csr_blocks );
844
845 for( i = 0; i < num_blocks; i++ )
846 {
847 hypre_TFree( hypre_ParCSRMatrixOffd( blocks[i] ) );
848 hypre_ParCSRMatrixOffd( blocks[i] ) = csr_blocks[i];
849 }
850
851 hypre_TFree( csr_blocks );
852 hypre_TFree( col_block_num );
853 hypre_TFree( row_block_num );
854
855 /* update block col-maps */
856 for( int bi = 0; bi < nr; bi++ )
857 {
858 for( int bj = 0; bj < nc; bj++ )
859 {
860 hypre_ParCSRMatrix* block = blocks[bi * nc + bj];
861 hypre_CSRMatrix* block_offd = hypre_ParCSRMatrixOffd( block );
862 HYPRE_Int block_offd_cols = hypre_CSRMatrixNumCols( block_offd );
863
864 HYPRE_Int* block_col_map = hypre_TAlloc( HYPRE_Int, block_offd_cols );
865 for( i = j = 0; i < offd_cols; i++ )
866 {
867 HYPRE_Int bn = offd_col_block_num[i];
868 if( bn == bj )
869 {
870 block_col_map[j++] = offd_global_col[i];
871 }
872 }
873 hypre_assert( j == block_offd_cols );
874
875 hypre_ParCSRMatrixColMapOffd( block ) = block_col_map;
876 }
877 }
878
879 hypre_TFree( offd_global_col );
880 hypre_TFree( offd_col_block_num );
881
882 /* finish the new matrices, make them own all the stuff */
883 for( i = 0; i < num_blocks; i++ )
884 {
885 hypre_ParCSRMatrixSetNumNonzeros( blocks[i] );
886 hypre_MatvecCommPkgCreate( blocks[i] );
887
888 hypre_ParCSRMatrixOwnsData( blocks[i] ) = 1;
889
890 /* only the first block will own the row/col_starts */
891 hypre_ParCSRMatrixOwnsRowStarts( blocks[i] ) = !i;
892 hypre_ParCSRMatrixOwnsColStarts( blocks[i] ) = !i;
893 }
894 }
895
896 /* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
897 void hypre_CSRMatrixBooleanMatvec( hypre_CSRMatrix* A,
898 HYPRE_Bool alpha,
899 HYPRE_Bool* x,
900 HYPRE_Bool beta,
901 HYPRE_Bool* y )
902 {
903 /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
904 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
905 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
906 HYPRE_Int num_rows = hypre_CSRMatrixNumRows( A );
907
908 HYPRE_Int* A_rownnz = hypre_CSRMatrixRownnz( A );
909 HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz( A );
910
911 HYPRE_Bool* x_data = x;
912 HYPRE_Bool* y_data = y;
913
914 HYPRE_Bool temp, tempx;
915
916 HYPRE_Int i, jj;
917
918 HYPRE_Int m;
919
920 HYPRE_Real xpar = 0.7;
921
922 /*-----------------------------------------------------------------------
923 * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
924 *-----------------------------------------------------------------------*/
925
926 if( alpha == 0 )
927 {
928 #ifdef HYPRE_USING_OPENMP
929 #pragma omp parallel for private( i ) HYPRE_SMP_SCHEDULE
930 #endif
931 for( i = 0; i < num_rows; i++ )
932 {
933 y_data[i] = y_data[i] && beta;
934 }
935 return;
936 }
937
938 /*-----------------------------------------------------------------------
939 * y = (beta/alpha)*y
940 *-----------------------------------------------------------------------*/
941
942 if( beta == 0 )
943 {
944 #ifdef HYPRE_USING_OPENMP
945 #pragma omp parallel for private( i ) HYPRE_SMP_SCHEDULE
946 #endif
947 for( i = 0; i < num_rows; i++ )
948 {
949 y_data[i] = 0;
950 }
951 }
952 else
953 {
954 /* beta is true -> no change to y_data */
955 }
956
957 /*-----------------------------------------------------------------
958 * y += A*x
959 *-----------------------------------------------------------------*/
960
961 /* use rownnz pointer to do the A*x multiplication when num_rownnz is smaller than num_rows
962 */
963
964 if( num_rownnz < xpar * ( num_rows ) )
965 {
966 #ifdef HYPRE_USING_OPENMP
967 #pragma omp parallel for private( i, jj, m, tempx ) HYPRE_SMP_SCHEDULE
968 #endif
969 for( i = 0; i < num_rownnz; i++ )
970 {
971 m = A_rownnz[i];
972
973 tempx = 0;
974 for( jj = A_i[m]; jj < A_i[m + 1]; jj++ )
975 {
976 /* tempx = tempx || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
977 tempx = tempx || x_data[A_j[jj]];
978 }
979 y_data[m] = y_data[m] || tempx;
980 }
981 }
982 else
983 {
984 #ifdef HYPRE_USING_OPENMP
985 #pragma omp parallel for private( i, jj, temp ) HYPRE_SMP_SCHEDULE
986 #endif
987 for( i = 0; i < num_rows; i++ )
988 {
989 temp = 0;
990 for( jj = A_i[i]; jj < A_i[i + 1]; jj++ )
991 {
992 /* temp = temp || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
993 temp = temp || x_data[A_j[jj]];
994 }
995 y_data[i] = y_data[i] || temp;
996 }
997 }
998
999 /*-----------------------------------------------------------------
1000 * y = alpha*y
1001 *-----------------------------------------------------------------*/
1002 /* alpha is true */
1003 }
1004
1005 /* Based on hypre_ParCSRCommHandleCreate in hypre's par_csr_communication.c. The
1006 input variable job controls the communication type: 1=Matvec, 2=MatvecT. */
1007 hypre_ParCSRCommHandle* hypre_ParCSRCommHandleCreate_bool( HYPRE_Int job,
1008 hypre_ParCSRCommPkg* comm_pkg,
1009 HYPRE_Bool* send_data,
1010 HYPRE_Bool* recv_data )
1011 {
1012 HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
1013 HYPRE_Int num_recvs = hypre_ParCSRCommPkgNumRecvs( comm_pkg );
1014 MPI_Comm comm = hypre_ParCSRCommPkgComm( comm_pkg );
1015
1016 hypre_ParCSRCommHandle* comm_handle;
1017 HYPRE_Int num_requests;
1018 hypre_MPI_Request* requests;
1019
1020 HYPRE_Int i, j;
1021 HYPRE_Int my_id, num_procs;
1022 HYPRE_Int ip, vec_start, vec_len;
1023
1024 num_requests = num_sends + num_recvs;
1025 requests = hypre_CTAlloc( hypre_MPI_Request, num_requests );
1026
1027 hypre_MPI_Comm_size( comm, &num_procs );
1028 hypre_MPI_Comm_rank( comm, &my_id );
1029
1030 j = 0;
1031 switch( job )
1032 {
1033 case 1: {
1034 HYPRE_Bool* d_send_data = (HYPRE_Bool*)send_data;
1035 HYPRE_Bool* d_recv_data = (HYPRE_Bool*)recv_data;
1036 for( i = 0; i < num_recvs; i++ )
1037 {
1038 ip = hypre_ParCSRCommPkgRecvProc( comm_pkg, i );
1039 vec_start = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i );
1040 vec_len = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i + 1 ) - vec_start;
1041 hypre_MPI_Irecv( &d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
1042 }
1043 for( i = 0; i < num_sends; i++ )
1044 {
1045 vec_start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
1046 vec_len = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ) - vec_start;
1047 ip = hypre_ParCSRCommPkgSendProc( comm_pkg, i );
1048 hypre_MPI_Isend( &d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
1049 }
1050 break;
1051 }
1052 case 2: {
1053 HYPRE_Bool* d_send_data = (HYPRE_Bool*)send_data;
1054 HYPRE_Bool* d_recv_data = (HYPRE_Bool*)recv_data;
1055 for( i = 0; i < num_sends; i++ )
1056 {
1057 vec_start = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
1058 vec_len = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ) - vec_start;
1059 ip = hypre_ParCSRCommPkgSendProc( comm_pkg, i );
1060 hypre_MPI_Irecv( &d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
1061 }
1062 for( i = 0; i < num_recvs; i++ )
1063 {
1064 ip = hypre_ParCSRCommPkgRecvProc( comm_pkg, i );
1065 vec_start = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i );
1066 vec_len = hypre_ParCSRCommPkgRecvVecStart( comm_pkg, i + 1 ) - vec_start;
1067 hypre_MPI_Isend( &d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL, ip, 0, comm, &requests[j++] );
1068 }
1069 break;
1070 }
1071 }
1072 /*--------------------------------------------------------------------
1073 * set up comm_handle and return
1074 *--------------------------------------------------------------------*/
1075
1076 comm_handle = hypre_CTAlloc( hypre_ParCSRCommHandle, 1 );
1077
1078 hypre_ParCSRCommHandleCommPkg( comm_handle ) = comm_pkg;
1079 hypre_ParCSRCommHandleSendData( comm_handle ) = send_data;
1080 hypre_ParCSRCommHandleRecvData( comm_handle ) = recv_data;
1081 hypre_ParCSRCommHandleNumRequests( comm_handle ) = num_requests;
1082 hypre_ParCSRCommHandleRequests( comm_handle ) = requests;
1083
1084 return comm_handle;
1085 }
1086
1087 /* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
1088 void hypre_ParCSRMatrixBooleanMatvec( hypre_ParCSRMatrix* A,
1089 HYPRE_Bool alpha,
1090 HYPRE_Bool* x,
1091 HYPRE_Bool beta,
1092 HYPRE_Bool* y )
1093 {
1094 hypre_ParCSRCommHandle* comm_handle;
1095 hypre_ParCSRCommPkg* comm_pkg = hypre_ParCSRMatrixCommPkg( A );
1096 hypre_CSRMatrix* diag = hypre_ParCSRMatrixDiag( A );
1097 hypre_CSRMatrix* offd = hypre_ParCSRMatrixOffd( A );
1098
1099 HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols( offd );
1100 HYPRE_Int num_sends, i, j, index;
1101
1102 HYPRE_Bool *x_tmp, *x_buf;
1103
1104 x_tmp = hypre_CTAlloc( HYPRE_Bool, num_cols_offd );
1105
1106 /*---------------------------------------------------------------------
1107 * If there exists no CommPkg for A, a CommPkg is generated using
1108 * equally load balanced partitionings
1109 *--------------------------------------------------------------------*/
1110 if( !comm_pkg )
1111 {
1112 hypre_MatvecCommPkgCreate( A );
1113 comm_pkg = hypre_ParCSRMatrixCommPkg( A );
1114 }
1115
1116 num_sends = hypre_ParCSRCommPkgNumSends( comm_pkg );
1117 x_buf = hypre_CTAlloc( HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart( comm_pkg, num_sends ) );
1118
1119 index = 0;
1120 for( i = 0; i < num_sends; i++ )
1121 {
1122 j = hypre_ParCSRCommPkgSendMapStart( comm_pkg, i );
1123 for( ; j < hypre_ParCSRCommPkgSendMapStart( comm_pkg, i + 1 ); j++ )
1124 {
1125 x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt( comm_pkg, j )];
1126 }
1127 }
1128
1129 comm_handle = hypre_ParCSRCommHandleCreate_bool( 1, comm_pkg, x_buf, x_tmp );
1130
1131 hypre_CSRMatrixBooleanMatvec( diag, alpha, x, beta, y );
1132
1133 hypre_ParCSRCommHandleDestroy( comm_handle );
1134
1135 if( num_cols_offd )
1136 {
1137 hypre_CSRMatrixBooleanMatvec( offd, alpha, x_tmp, 1, y );
1138 }
1139
1140 hypre_TFree( x_buf );
1141 hypre_TFree( x_tmp );
1142 }
1143
1144 HYPRE_Int hypre_CSRMatrixSum( hypre_CSRMatrix* A, HYPRE_Complex beta, hypre_CSRMatrix* B )
1145 {
1146 HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
1147 HYPRE_Int* A_i = hypre_CSRMatrixI( A );
1148 HYPRE_Int* A_j = hypre_CSRMatrixJ( A );
1149 HYPRE_Int nrows_A = hypre_CSRMatrixNumRows( A );
1150 HYPRE_Int ncols_A = hypre_CSRMatrixNumCols( A );
1151 HYPRE_Complex* B_data = hypre_CSRMatrixData( B );
1152 HYPRE_Int* B_i = hypre_CSRMatrixI( B );
1153 HYPRE_Int* B_j = hypre_CSRMatrixJ( B );
1154 HYPRE_Int nrows_B = hypre_CSRMatrixNumRows( B );
1155 HYPRE_Int ncols_B = hypre_CSRMatrixNumCols( B );
1156
1157 HYPRE_Int ia, j, pos;
1158 HYPRE_Int* marker;
1159
1160 if( nrows_A != nrows_B || ncols_A != ncols_B )
1161 {
1162 return -1; /* error: incompatible matrix dimensions */
1163 }
1164
1165 marker = hypre_CTAlloc( HYPRE_Int, ncols_A );
1166 for( ia = 0; ia < ncols_A; ia++ )
1167 {
1168 marker[ia] = -1;
1169 }
1170
1171 for( ia = 0; ia < nrows_A; ia++ )
1172 {
1173 for( j = A_i[ia]; j < A_i[ia + 1]; j++ )
1174 {
1175 marker[A_j[j]] = j;
1176 }
1177
1178 for( j = B_i[ia]; j < B_i[ia + 1]; j++ )
1179 {
1180 pos = marker[B_j[j]];
1181 if( pos < A_i[ia] )
1182 {
1183 return -2; /* error: found an entry in B that is not present in A */
1184 }
1185 A_data[pos] += beta * B_data[j];
1186 }
1187 }
1188
1189 hypre_TFree( marker );
1190 return 0;
1191 }
1192
1193 hypre_ParCSRMatrix* hypre_ParCSRMatrixAdd( hypre_ParCSRMatrix* A, hypre_ParCSRMatrix* B )
1194 {
1195 MPI_Comm comm = hypre_ParCSRMatrixComm( A );
1196 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
1197 hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
1198 HYPRE_Int* A_cmap = hypre_ParCSRMatrixColMapOffd( A );
1199 HYPRE_Int A_cmap_size = hypre_CSRMatrixNumCols( A_offd );
1200 hypre_CSRMatrix* B_diag = hypre_ParCSRMatrixDiag( B );
1201 hypre_CSRMatrix* B_offd = hypre_ParCSRMatrixOffd( B );
1202 HYPRE_Int* B_cmap = hypre_ParCSRMatrixColMapOffd( B );
1203 HYPRE_Int B_cmap_size = hypre_CSRMatrixNumCols( B_offd );
1204 hypre_ParCSRMatrix* C;
1205 hypre_CSRMatrix* C_diag;
1206 hypre_CSRMatrix* C_offd;
1207 HYPRE_Int* C_cmap;
1208 HYPRE_Int im;
1209 HYPRE_Int cmap_differ;
1210
1211 /* Check if A_cmap and B_cmap are the same. */
1212 cmap_differ = 0;
1213 if( A_cmap_size != B_cmap_size )
1214 {
1215 cmap_differ = 1; /* A and B have different cmap_size */
1216 }
1217 else
1218 {
1219 for( im = 0; im < A_cmap_size; im++ )
1220 {
1221 if( A_cmap[im] != B_cmap[im] )
1222 {
1223 cmap_differ = 1; /* A and B have different cmap arrays */
1224 break;
1225 }
1226 }
1227 }
1228
1229 if( cmap_differ == 0 )
1230 {
1231 /* A and B have the same column mapping for their off-diagonal blocks so
1232 we can sum the diagonal and off-diagonal blocks separately and reduce
1233 temporary memory usage. */
1234
1235 /* Add diagonals, off-diagonals, copy cmap. */
1236 C_diag = hypre_CSRMatrixAdd( A_diag, B_diag );
1237 if( !C_diag )
1238 {
1239 return NULL; /* error: A_diag and B_diag have different dimensions */
1240 }
1241 C_offd = hypre_CSRMatrixAdd( A_offd, B_offd );
1242 if( !C_offd )
1243 {
1244 hypre_CSRMatrixDestroy( C_diag );
1245 return NULL; /* error: A_offd and B_offd have different dimensions */
1246 }
1247 /* copy A_cmap -> C_cmap */
1248 C_cmap = hypre_TAlloc( HYPRE_Int, A_cmap_size );
1249 for( im = 0; im < A_cmap_size; im++ )
1250 {
1251 C_cmap[im] = A_cmap[im];
1252 }
1253
1254 C = hypre_ParCSRMatrixCreate( comm, hypre_ParCSRMatrixGlobalNumRows( A ),
1255 hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
1256 hypre_ParCSRMatrixColStarts( A ), hypre_CSRMatrixNumCols( C_offd ),
1257 hypre_CSRMatrixNumNonzeros( C_diag ), hypre_CSRMatrixNumNonzeros( C_offd ) );
1258
1259 /* In C, destroy diag/offd (allocated by Create) and replace them with
1260 C_diag/C_offd. */
1261 hypre_CSRMatrixDestroy( hypre_ParCSRMatrixDiag( C ) );
1262 hypre_CSRMatrixDestroy( hypre_ParCSRMatrixOffd( C ) );
1263 hypre_ParCSRMatrixDiag( C ) = C_diag;
1264 hypre_ParCSRMatrixOffd( C ) = C_offd;
1265
1266 hypre_ParCSRMatrixColMapOffd( C ) = C_cmap;
1267 }
1268 else
1269 {
1270 /* A and B have different column mappings for their off-diagonal blocks so
1271 we need to use the column maps to create full-width CSR matricies. */
1272
1273 int ierr = 0;
1274 hypre_CSRMatrix* csr_A;
1275 hypre_CSRMatrix* csr_B;
1276 hypre_CSRMatrix* csr_C_temp;
1277
1278 /* merge diag and off-diag portions of A */
1279 csr_A = hypre_MergeDiagAndOffd( A );
1280
1281 /* merge diag and off-diag portions of B */
1282 csr_B = hypre_MergeDiagAndOffd( B );
1283
1284 /* add A and B */
1285 csr_C_temp = hypre_CSRMatrixAdd( csr_A, csr_B );
1286
1287 /* delete CSR versions of A and B */
1288 ierr += hypre_CSRMatrixDestroy( csr_A );
1289 ierr += hypre_CSRMatrixDestroy( csr_B );
1290
1291 /* create a new empty ParCSR matrix to contain the sum */
1292 C = hypre_ParCSRMatrixCreate( hypre_ParCSRMatrixComm( A ), hypre_ParCSRMatrixGlobalNumRows( A ),
1293 hypre_ParCSRMatrixGlobalNumCols( A ), hypre_ParCSRMatrixRowStarts( A ),
1294 hypre_ParCSRMatrixColStarts( A ), 0, 0, 0 );
1295
1296 /* split C into diag and off-diag portions */
1297 /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the
1298 number of columns in csr_C_temp which is the global number of columns
1299 in A and B. This does not scale well. */
1300 ierr += GenerateDiagAndOffd( csr_C_temp, C, hypre_ParCSRMatrixFirstColDiag( A ),
1301 hypre_ParCSRMatrixLastColDiag( A ) );
1302
1303 /* delete CSR version of C */
1304 ierr += hypre_CSRMatrixDestroy( csr_C_temp );
1305 }
1306
1307 /* hypre_ParCSRMatrixSetNumNonzeros(A); */
1308
1309 /* Make sure that the first entry in each row is the diagonal one. */
1310 hypre_CSRMatrixReorder( hypre_ParCSRMatrixDiag( C ) );
1311
1312 /* C owns diag, offd, and cmap. */
1313 hypre_ParCSRMatrixSetDataOwner( C, 1 );
1314 /* C does not own row and column starts. */
1315 hypre_ParCSRMatrixSetRowStartsOwner( C, 0 );
1316 hypre_ParCSRMatrixSetColStartsOwner( C, 0 );
1317
1318 return C;
1319 }
1320
1321 HYPRE_Int hypre_ParCSRMatrixSum( hypre_ParCSRMatrix* A, HYPRE_Complex beta, hypre_ParCSRMatrix* B )
1322 {
1323 hypre_CSRMatrix* A_diag = hypre_ParCSRMatrixDiag( A );
1324 hypre_CSRMatrix* A_offd = hypre_ParCSRMatrixOffd( A );
1325 hypre_CSRMatrix* B_diag = hypre_ParCSRMatrixDiag( B );
1326 hypre_CSRMatrix* B_offd = hypre_ParCSRMatrixOffd( B );
1327 HYPRE_Int error;
1328
1329 error = hypre_CSRMatrixSum( A_diag, beta, B_diag );
1330 error = error ? error : hypre_CSRMatrixSum( A_offd, beta, B_offd );
1331
1332 return error;
1333 }
1334
1335 HYPRE_Int hypre_CSRMatrixSetConstantValues( hypre_CSRMatrix* A, HYPRE_Complex value )
1336 {
1337 HYPRE_Complex* A_data = hypre_CSRMatrixData( A );
1338 HYPRE_Int A_nnz = hypre_CSRMatrixNumNonzeros( A );
1339 HYPRE_Int ia;
1340
1341 for( ia = 0; ia < A_nnz; ia++ )
1342 {
1343 A_data[ia] = value;
1344 }
1345
1346 return 0;
1347 }
1348
1349 HYPRE_Int hypre_ParCSRMatrixSetConstantValues( hypre_ParCSRMatrix* A, HYPRE_Complex value )
1350 {
1351 hypre_CSRMatrixSetConstantValues( hypre_ParCSRMatrixDiag( A ), value );
1352 hypre_CSRMatrixSetConstantValues( hypre_ParCSRMatrixOffd( A ), value );
1353
1354 return 0;
1355 }
1356
1357 } // namespace internal
1358
1359 } // namespace moab
1360
1361 #endif // MOAB_HAVE_MPI
1.9.1.