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