Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
hypre_parcsr.cpp
Go to the documentation of this file.
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