MFEM  v4.1.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-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
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  Eliminate rows of A, setting all entries in the eliminated rows to zero.
464 */
465 void hypre_CSRMatrixEliminateRows(hypre_CSRMatrix *A,
466  HYPRE_Int nrows, const HYPRE_Int *rows)
467 {
468  HYPRE_Int irow, i, j;
469  HYPRE_Int A_beg, A_end;
470 
471  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
472  HYPRE_Real *A_data = hypre_CSRMatrixData(A);
473 
474  for (i = 0; i < nrows; i++)
475  {
476  irow = rows[i];
477  A_beg = A_i[irow];
478  A_end = A_i[irow+1];
479  /* eliminate row */
480  for (j = A_beg; j < A_end; j++)
481  {
482  A_data[j] = 0.0;
483  }
484  }
485 }
486 
487 
488 /*
489  Function: hypre_ParCSRMatrixEliminateAAe
490 
491  (input) (output)
492 
493  / A_ii | A_ib \ / A_ii | 0 \
494  A = | -----+----- | ---> | -----+----- |
495  \ A_bi | A_bb / \ 0 | I /
496 
497 
498  / 0 | A_ib \
499  Ae = | -----+--------- |
500  \ A_bi | A_bb - I /
501 
502 */
503 void hypre_ParCSRMatrixEliminateAAe(hypre_ParCSRMatrix *A,
504  hypre_ParCSRMatrix **Ae,
505  HYPRE_Int num_rowscols_to_elim,
506  HYPRE_Int *rowscols_to_elim,
507  int ignore_rows)
508 {
509  HYPRE_Int i, j, k;
510 
511  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
512  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
513  HYPRE_Int A_diag_ncols = hypre_CSRMatrixNumCols(A_diag);
514  HYPRE_Int A_offd_ncols = hypre_CSRMatrixNumCols(A_offd);
515 
516  *Ae = hypre_ParCSRMatrixCreate(hypre_ParCSRMatrixComm(A),
517  hypre_ParCSRMatrixGlobalNumRows(A),
518  hypre_ParCSRMatrixGlobalNumCols(A),
519  hypre_ParCSRMatrixRowStarts(A),
520  hypre_ParCSRMatrixColStarts(A),
521  0, 0, 0);
522 
523  hypre_ParCSRMatrixSetRowStartsOwner(*Ae, 0);
524  hypre_ParCSRMatrixSetColStartsOwner(*Ae, 0);
525 
526  hypre_CSRMatrix *Ae_diag = hypre_ParCSRMatrixDiag(*Ae);
527  hypre_CSRMatrix *Ae_offd = hypre_ParCSRMatrixOffd(*Ae);
528  HYPRE_Int Ae_offd_ncols;
529 
530  HYPRE_Int num_offd_cols_to_elim;
531  HYPRE_Int *offd_cols_to_elim;
532 
533  HYPRE_Int *A_col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
534  HYPRE_Int *Ae_col_map_offd;
535 
536  HYPRE_Int *col_mark;
537  HYPRE_Int *col_remap;
538 
539  /* figure out which offd cols should be eliminated */
540  {
541  hypre_ParCSRCommHandle *comm_handle;
542  hypre_ParCSRCommPkg *comm_pkg;
543  HYPRE_Int num_sends, *int_buf_data;
544  HYPRE_Int index, start;
545 
546  HYPRE_Int *eliminate_diag_col = mfem_hypre_CTAlloc(HYPRE_Int, A_diag_ncols);
547  HYPRE_Int *eliminate_offd_col = mfem_hypre_CTAlloc(HYPRE_Int, A_offd_ncols);
548 
549  /* make sure A has a communication package */
550  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
551  if (!comm_pkg)
552  {
553  hypre_MatvecCommPkgCreate(A);
554  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
555  }
556 
557  /* which of the local rows are to be eliminated */
558  for (i = 0; i < A_diag_ncols; i++)
559  {
560  eliminate_diag_col[i] = 0;
561  }
562  for (i = 0; i < num_rowscols_to_elim; i++)
563  {
564  eliminate_diag_col[rowscols_to_elim[i]] = 1;
565  }
566 
567  /* use a Matvec communication pattern to find (in eliminate_col)
568  which of the local offd columns are to be eliminated */
569  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
570  int_buf_data = mfem_hypre_CTAlloc(
571  HYPRE_Int,
572  hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
573  index = 0;
574  for (i = 0; i < num_sends; i++)
575  {
576  start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
577  for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
578  {
579  k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
580  int_buf_data[index++] = eliminate_diag_col[k];
581  }
582  }
583  comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg,
584  int_buf_data, eliminate_offd_col);
585 
586  /* eliminate diagonal part, overlapping it with communication */
587  if (ignore_rows)
588  {
589  hypre_CSRMatrixElimCreate(A_diag, Ae_diag,
590  0, nullptr,
591  num_rowscols_to_elim, rowscols_to_elim,
592  NULL);
593 
594  hypre_CSRMatrixEliminateRowsCols(A_diag, Ae_diag,
595  0, nullptr,
596  num_rowscols_to_elim, rowscols_to_elim,
597  1, NULL);
598  }
599  else
600  {
601  hypre_CSRMatrixElimCreate(A_diag, Ae_diag,
602  num_rowscols_to_elim, rowscols_to_elim,
603  num_rowscols_to_elim, rowscols_to_elim,
604  NULL);
605 
606  hypre_CSRMatrixEliminateRowsCols(A_diag, Ae_diag,
607  num_rowscols_to_elim, rowscols_to_elim,
608  num_rowscols_to_elim, rowscols_to_elim,
609  1, NULL);
610  }
611 
612  hypre_CSRMatrixReorder(Ae_diag);
613 
614  /* finish the communication */
615  hypre_ParCSRCommHandleDestroy(comm_handle);
616 
617  /* received eliminate_col[], count offd columns to eliminate */
618  num_offd_cols_to_elim = 0;
619  for (i = 0; i < A_offd_ncols; i++)
620  {
621  if (eliminate_offd_col[i]) { num_offd_cols_to_elim++; }
622  }
623 
624  offd_cols_to_elim = mfem_hypre_CTAlloc(HYPRE_Int, num_offd_cols_to_elim);
625 
626  /* get a list of offd column indices and coefs */
627  num_offd_cols_to_elim = 0;
628  for (i = 0; i < A_offd_ncols; i++)
629  {
630  if (eliminate_offd_col[i])
631  {
632  offd_cols_to_elim[num_offd_cols_to_elim++] = i;
633  }
634  }
635 
636  mfem_hypre_TFree(int_buf_data);
637  mfem_hypre_TFree(eliminate_offd_col);
638  mfem_hypre_TFree(eliminate_diag_col);
639  }
640 
641  /* eliminate the off-diagonal part */
642  col_mark = mfem_hypre_CTAlloc(HYPRE_Int, A_offd_ncols);
643  col_remap = mfem_hypre_CTAlloc(HYPRE_Int, A_offd_ncols);
644 
645  if (ignore_rows)
646  {
647  hypre_CSRMatrixElimCreate(A_offd, Ae_offd,
648  0, nullptr,
649  num_offd_cols_to_elim, offd_cols_to_elim,
650  col_mark);
651 
652  for (i = k = 0; i < A_offd_ncols; i++)
653  {
654  if (col_mark[i]) { col_remap[i] = k++; }
655  }
656 
657  hypre_CSRMatrixEliminateRowsCols(A_offd, Ae_offd,
658  0, nullptr,
659  num_offd_cols_to_elim, offd_cols_to_elim,
660  0, col_remap);
661  }
662  else
663  {
664  hypre_CSRMatrixElimCreate(A_offd, Ae_offd,
665  num_rowscols_to_elim, rowscols_to_elim,
666  num_offd_cols_to_elim, offd_cols_to_elim,
667  col_mark);
668 
669  for (i = k = 0; i < A_offd_ncols; i++)
670  {
671  if (col_mark[i]) { col_remap[i] = k++; }
672  }
673 
674  hypre_CSRMatrixEliminateRowsCols(A_offd, Ae_offd,
675  num_rowscols_to_elim, rowscols_to_elim,
676  num_offd_cols_to_elim, offd_cols_to_elim,
677  0, col_remap);
678  }
679 
680  /* create col_map_offd for Ae */
681  Ae_offd_ncols = 0;
682  for (i = 0; i < A_offd_ncols; i++)
683  {
684  if (col_mark[i]) { Ae_offd_ncols++; }
685  }
686 
687  Ae_col_map_offd = mfem_hypre_CTAlloc(HYPRE_Int, Ae_offd_ncols);
688 
689  Ae_offd_ncols = 0;
690  for (i = 0; i < A_offd_ncols; i++)
691  {
692  if (col_mark[i])
693  {
694  Ae_col_map_offd[Ae_offd_ncols++] = A_col_map_offd[i];
695  }
696  }
697 
698  hypre_ParCSRMatrixColMapOffd(*Ae) = Ae_col_map_offd;
699  hypre_CSRMatrixNumCols(Ae_offd) = Ae_offd_ncols;
700 
701  mfem_hypre_TFree(col_remap);
702  mfem_hypre_TFree(col_mark);
703  mfem_hypre_TFree(offd_cols_to_elim);
704 
705  hypre_ParCSRMatrixSetNumNonzeros(*Ae);
706  hypre_MatvecCommPkgCreate(*Ae);
707 }
708 
709 
710 // Eliminate rows from the diagonal and off-diagonal blocks of the matrix
711 void hypre_ParCSRMatrixEliminateRows(hypre_ParCSRMatrix *A,
712  HYPRE_Int num_rows_to_elim,
713  const HYPRE_Int *rows_to_elim)
714 {
715  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
716  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
717  hypre_CSRMatrixEliminateRows(A_diag, num_rows_to_elim, rows_to_elim);
718  hypre_CSRMatrixEliminateRows(A_offd, num_rows_to_elim, rows_to_elim);
719 }
720 
721 
722 /*--------------------------------------------------------------------------
723  * Split
724  *--------------------------------------------------------------------------*/
725 
726 void hypre_CSRMatrixSplit(hypre_CSRMatrix *A,
727  HYPRE_Int nr, HYPRE_Int nc,
728  HYPRE_Int *row_block_num, HYPRE_Int *col_block_num,
729  hypre_CSRMatrix **blocks)
730 {
731  HYPRE_Int i, j, k, bi, bj;
732 
733  HYPRE_Int* A_i = hypre_CSRMatrixI(A);
734  HYPRE_Int* A_j = hypre_CSRMatrixJ(A);
735  HYPRE_Complex* A_data = hypre_CSRMatrixData(A);
736 
737  HYPRE_Int A_rows = hypre_CSRMatrixNumRows(A);
738  HYPRE_Int A_cols = hypre_CSRMatrixNumCols(A);
739 
740  HYPRE_Int *num_rows = mfem_hypre_CTAlloc(HYPRE_Int, nr);
741  HYPRE_Int *num_cols = mfem_hypre_CTAlloc(HYPRE_Int, nc);
742 
743  HYPRE_Int *block_row = mfem_hypre_TAlloc(HYPRE_Int, A_rows);
744  HYPRE_Int *block_col = mfem_hypre_TAlloc(HYPRE_Int, A_cols);
745 
746  for (i = 0; i < A_rows; i++)
747  {
748  block_row[i] = num_rows[row_block_num[i]]++;
749  }
750  for (j = 0; j < A_cols; j++)
751  {
752  block_col[j] = num_cols[col_block_num[j]]++;
753  }
754 
755  /* allocate the blocks */
756  for (i = 0; i < nr; i++)
757  {
758  for (j = 0; j < nc; j++)
759  {
760  hypre_CSRMatrix *B = hypre_CSRMatrixCreate(num_rows[i], num_cols[j], 0);
761  hypre_CSRMatrixI(B) = mfem_hypre_CTAlloc(HYPRE_Int, num_rows[i] + 1);
762  blocks[i*nc + j] = B;
763  }
764  }
765 
766  /* count block row nnz */
767  for (i = 0; i < A_rows; i++)
768  {
769  bi = row_block_num[i];
770  for (j = A_i[i]; j < A_i[i+1]; j++)
771  {
772  bj = col_block_num[A_j[j]];
773  hypre_CSRMatrix *B = blocks[bi*nc + bj];
774  hypre_CSRMatrixI(B)[block_row[i] + 1]++;
775  }
776  }
777 
778  /* count block nnz */
779  for (k = 0; k < nr*nc; k++)
780  {
781  hypre_CSRMatrix *B = blocks[k];
782  HYPRE_Int* B_i = hypre_CSRMatrixI(B);
783 
784  HYPRE_Int nnz = 0, rs;
785  for (int k = 1; k <= hypre_CSRMatrixNumRows(B); k++)
786  {
787  rs = B_i[k], B_i[k] = nnz, nnz += rs;
788  }
789 
790  hypre_CSRMatrixJ(B) = mfem_hypre_TAlloc(HYPRE_Int, nnz);
791  hypre_CSRMatrixData(B) = mfem_hypre_TAlloc(HYPRE_Complex, nnz);
792  hypre_CSRMatrixNumNonzeros(B) = nnz;
793  }
794 
795  /* populate blocks */
796  for (i = 0; i < A_rows; i++)
797  {
798  bi = row_block_num[i];
799  for (j = A_i[i]; j < A_i[i+1]; j++)
800  {
801  k = A_j[j];
802  bj = col_block_num[k];
803  hypre_CSRMatrix *B = blocks[bi*nc + bj];
804  HYPRE_Int *bii = hypre_CSRMatrixI(B) + block_row[i] + 1;
805  hypre_CSRMatrixJ(B)[*bii] = block_col[k];
806  hypre_CSRMatrixData(B)[*bii] = A_data[j];
807  (*bii)++;
808  }
809  }
810 
811  mfem_hypre_TFree(block_col);
812  mfem_hypre_TFree(block_row);
813 
814  mfem_hypre_TFree(num_cols);
815  mfem_hypre_TFree(num_rows);
816 }
817 
818 
819 void hypre_ParCSRMatrixSplit(hypre_ParCSRMatrix *A,
820  HYPRE_Int nr, HYPRE_Int nc,
821  hypre_ParCSRMatrix **blocks,
822  int interleaved_rows, int interleaved_cols)
823 {
824  HYPRE_Int i, j, k;
825 
826  MPI_Comm comm = hypre_ParCSRMatrixComm(A);
827 
828  hypre_CSRMatrix *Adiag = hypre_ParCSRMatrixDiag(A);
829  hypre_CSRMatrix *Aoffd = hypre_ParCSRMatrixOffd(A);
830 
831  HYPRE_Int global_rows = hypre_ParCSRMatrixGlobalNumRows(A);
832  HYPRE_Int global_cols = hypre_ParCSRMatrixGlobalNumCols(A);
833 
834  HYPRE_Int local_rows = hypre_CSRMatrixNumRows(Adiag);
835  HYPRE_Int local_cols = hypre_CSRMatrixNumCols(Adiag);
836  HYPRE_Int offd_cols = hypre_CSRMatrixNumCols(Aoffd);
837 
838  hypre_assert(local_rows % nr == 0 && local_cols % nc == 0);
839  hypre_assert(global_rows % nr == 0 && global_cols % nc == 0);
840 
841  HYPRE_Int block_rows = local_rows / nr;
842  HYPRE_Int block_cols = local_cols / nc;
843  HYPRE_Int num_blocks = nr * nc;
844 
845  /* mark local rows and columns with block number */
846  HYPRE_Int *row_block_num = mfem_hypre_TAlloc(HYPRE_Int, local_rows);
847  HYPRE_Int *col_block_num = mfem_hypre_TAlloc(HYPRE_Int, local_cols);
848 
849  for (i = 0; i < local_rows; i++)
850  {
851  row_block_num[i] = interleaved_rows ? (i % nr) : (i / block_rows);
852  }
853  for (i = 0; i < local_cols; i++)
854  {
855  col_block_num[i] = interleaved_cols ? (i % nc) : (i / block_cols);
856  }
857 
858  /* determine the block numbers for offd columns */
859  HYPRE_Int* offd_col_block_num = mfem_hypre_TAlloc(HYPRE_Int, offd_cols);
860  hypre_ParCSRCommHandle *comm_handle;
861  HYPRE_Int *int_buf_data;
862  {
863  /* make sure A has a communication package */
864  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
865  if (!comm_pkg)
866  {
867  hypre_MatvecCommPkgCreate(A);
868  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
869  }
870 
871  /* calculate the final global column numbers for each block */
872  HYPRE_Int *count = mfem_hypre_CTAlloc(HYPRE_Int, nc);
873  HYPRE_Int *block_global_col = mfem_hypre_TAlloc(HYPRE_Int, local_cols);
874  HYPRE_Int first_col = hypre_ParCSRMatrixFirstColDiag(A) / nc;
875  for (i = 0; i < local_cols; i++)
876  {
877  block_global_col[i] = first_col + count[col_block_num[i]]++;
878  }
879  mfem_hypre_TFree(count);
880 
881  /* use a Matvec communication pattern to determine offd_col_block_num */
882  HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
883  int_buf_data = mfem_hypre_CTAlloc(
884  HYPRE_Int,
885  hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
886  HYPRE_Int start, index = 0;
887  for (i = 0; i < num_sends; i++)
888  {
889  start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
890  for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
891  {
892  k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
893  int_buf_data[index++] = col_block_num[k] + nc*block_global_col[k];
894  }
895  }
896  mfem_hypre_TFree(block_global_col);
897 
898  comm_handle = hypre_ParCSRCommHandleCreate(11, comm_pkg, int_buf_data,
899  offd_col_block_num);
900  }
901 
902  /* create the block matrices */
903  HYPRE_Int num_procs = 1;
904  if (!hypre_ParCSRMatrixAssumedPartition(A))
905  {
906  hypre_MPI_Comm_size(comm, &num_procs);
907  }
908 
909  HYPRE_Int *row_starts = mfem_hypre_TAlloc(HYPRE_Int, num_procs+1);
910  HYPRE_Int *col_starts = mfem_hypre_TAlloc(HYPRE_Int, num_procs+1);
911  for (i = 0; i <= num_procs; i++)
912  {
913  row_starts[i] = hypre_ParCSRMatrixRowStarts(A)[i] / nr;
914  col_starts[i] = hypre_ParCSRMatrixColStarts(A)[i] / nc;
915  }
916 
917  for (i = 0; i < num_blocks; i++)
918  {
919  blocks[i] = hypre_ParCSRMatrixCreate(comm,
920  global_rows/nr, global_cols/nc,
921  row_starts, col_starts, 0, 0, 0);
922  }
923 
924  /* split diag part */
925  hypre_CSRMatrix **csr_blocks = mfem_hypre_TAlloc(hypre_CSRMatrix*, nr*nc);
926  hypre_CSRMatrixSplit(Adiag, nr, nc, row_block_num, col_block_num,
927  csr_blocks);
928 
929  for (i = 0; i < num_blocks; i++)
930  {
931  mfem_hypre_TFree(hypre_ParCSRMatrixDiag(blocks[i]));
932  hypre_ParCSRMatrixDiag(blocks[i]) = csr_blocks[i];
933  }
934 
935  /* finish communication, receive offd_col_block_num */
936  hypre_ParCSRCommHandleDestroy(comm_handle);
937  mfem_hypre_TFree(int_buf_data);
938 
939  /* decode global offd column numbers */
940  HYPRE_Int* offd_global_col = mfem_hypre_TAlloc(HYPRE_Int, offd_cols);
941  for (i = 0; i < offd_cols; i++)
942  {
943  offd_global_col[i] = offd_col_block_num[i] / nc;
944  offd_col_block_num[i] %= nc;
945  }
946 
947  /* split offd part */
948  hypre_CSRMatrixSplit(Aoffd, nr, nc, row_block_num, offd_col_block_num,
949  csr_blocks);
950 
951  for (i = 0; i < num_blocks; i++)
952  {
953  mfem_hypre_TFree(hypre_ParCSRMatrixOffd(blocks[i]));
954  hypre_ParCSRMatrixOffd(blocks[i]) = csr_blocks[i];
955  }
956 
957  mfem_hypre_TFree(csr_blocks);
958  mfem_hypre_TFree(col_block_num);
959  mfem_hypre_TFree(row_block_num);
960 
961  /* update block col-maps */
962  for (int bi = 0; bi < nr; bi++)
963  {
964  for (int bj = 0; bj < nc; bj++)
965  {
966  hypre_ParCSRMatrix *block = blocks[bi*nc + bj];
967  hypre_CSRMatrix *block_offd = hypre_ParCSRMatrixOffd(block);
968  HYPRE_Int block_offd_cols = hypre_CSRMatrixNumCols(block_offd);
969 
970  HYPRE_Int *block_col_map = mfem_hypre_TAlloc(HYPRE_Int,
971  block_offd_cols);
972  for (i = j = 0; i < offd_cols; i++)
973  {
974  HYPRE_Int bn = offd_col_block_num[i];
975  if (bn == bj) { block_col_map[j++] = offd_global_col[i]; }
976  }
977  hypre_assert(j == block_offd_cols);
978 
979  hypre_ParCSRMatrixColMapOffd(block) = block_col_map;
980  }
981  }
982 
983  mfem_hypre_TFree(offd_global_col);
984  mfem_hypre_TFree(offd_col_block_num);
985 
986  /* finish the new matrices, make them own all the stuff */
987  for (i = 0; i < num_blocks; i++)
988  {
989  hypre_ParCSRMatrixSetNumNonzeros(blocks[i]);
990  hypre_MatvecCommPkgCreate(blocks[i]);
991 
992  hypre_ParCSRMatrixOwnsData(blocks[i]) = 1;
993 
994  /* only the first block will own the row/col_starts */
995  hypre_ParCSRMatrixOwnsRowStarts(blocks[i]) = !i;
996  hypre_ParCSRMatrixOwnsColStarts(blocks[i]) = !i;
997  }
998 }
999 
1000 /* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
1001 void hypre_CSRMatrixBooleanMatvec(hypre_CSRMatrix *A,
1002  HYPRE_Bool alpha,
1003  HYPRE_Bool *x,
1004  HYPRE_Bool beta,
1005  HYPRE_Bool *y)
1006 {
1007  /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
1008  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1009  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1010  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1011 
1012  HYPRE_Int *A_rownnz = hypre_CSRMatrixRownnz(A);
1013  HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz(A);
1014 
1015  HYPRE_Bool *x_data = x;
1016  HYPRE_Bool *y_data = y;
1017 
1018  HYPRE_Bool temp, tempx;
1019 
1020  HYPRE_Int i, jj;
1021 
1022  HYPRE_Int m;
1023 
1024  HYPRE_Real xpar=0.7;
1025 
1026  /*-----------------------------------------------------------------------
1027  * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
1028  *-----------------------------------------------------------------------*/
1029 
1030  if (alpha == 0)
1031  {
1032 #ifdef HYPRE_USING_OPENMP
1033  #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
1034 #endif
1035  for (i = 0; i < num_rows; i++)
1036  {
1037  y_data[i] = y_data[i] && beta;
1038  }
1039  return;
1040  }
1041 
1042  /*-----------------------------------------------------------------------
1043  * y = (beta/alpha)*y
1044  *-----------------------------------------------------------------------*/
1045 
1046  if (beta == 0)
1047  {
1048 #ifdef HYPRE_USING_OPENMP
1049  #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
1050 #endif
1051  for (i = 0; i < num_rows; i++)
1052  {
1053  y_data[i] = 0;
1054  }
1055  }
1056  else
1057  {
1058  /* beta is true -> no change to y_data */
1059  }
1060 
1061  /*-----------------------------------------------------------------
1062  * y += A*x
1063  *-----------------------------------------------------------------*/
1064 
1065  /* use rownnz pointer to do the A*x multiplication when num_rownnz is smaller than num_rows */
1066 
1067  if (num_rownnz < xpar*(num_rows))
1068  {
1069 #ifdef HYPRE_USING_OPENMP
1070  #pragma omp parallel for private(i,jj,m,tempx) HYPRE_SMP_SCHEDULE
1071 #endif
1072  for (i = 0; i < num_rownnz; i++)
1073  {
1074  m = A_rownnz[i];
1075 
1076  tempx = 0;
1077  for (jj = A_i[m]; jj < A_i[m+1]; jj++)
1078  {
1079  /* tempx = tempx || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
1080  tempx = tempx || x_data[A_j[jj]];
1081  }
1082  y_data[m] = y_data[m] || tempx;
1083  }
1084  }
1085  else
1086  {
1087 #ifdef HYPRE_USING_OPENMP
1088  #pragma omp parallel for private(i,jj,temp) HYPRE_SMP_SCHEDULE
1089 #endif
1090  for (i = 0; i < num_rows; i++)
1091  {
1092  temp = 0;
1093  for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1094  {
1095  /* temp = temp || ((A_data[jj] != 0.0) && x_data[A_j[jj]]); */
1096  temp = temp || x_data[A_j[jj]];
1097  }
1098  y_data[i] = y_data[i] || temp;
1099  }
1100  }
1101 
1102  /*-----------------------------------------------------------------
1103  * y = alpha*y
1104  *-----------------------------------------------------------------*/
1105  /* alpha is true */
1106 }
1107 
1108 /* Based on hypre_CSRMatrixMatvecT in hypre's csr_matvec.c */
1109 void hypre_CSRMatrixBooleanMatvecT(hypre_CSRMatrix *A,
1110  HYPRE_Bool alpha,
1111  HYPRE_Bool *x,
1112  HYPRE_Bool beta,
1113  HYPRE_Bool *y)
1114 {
1115  /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
1116  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1117  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1118  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
1119  HYPRE_Int num_cols = hypre_CSRMatrixNumCols(A);
1120 
1121  HYPRE_Bool *x_data = x;
1122  HYPRE_Bool *y_data = y;
1123 
1124  HYPRE_Int i, j, jj;
1125 
1126  /*-----------------------------------------------------------------------
1127  * y = beta*y
1128  *-----------------------------------------------------------------------*/
1129 
1130  if (beta == 0)
1131  {
1132  for (i = 0; i < num_cols; i++)
1133  {
1134  y_data[i] = 0;
1135  }
1136  }
1137  else
1138  {
1139  /* beta is true -> no change to y_data */
1140  }
1141 
1142  /*-----------------------------------------------------------------------
1143  * Check if (alpha == 0)
1144  *-----------------------------------------------------------------------*/
1145 
1146  if (alpha == 0)
1147  {
1148  return;
1149  }
1150 
1151  /* alpha is true */
1152 
1153  /*-----------------------------------------------------------------
1154  * y += A^T*x
1155  *-----------------------------------------------------------------*/
1156  for (i = 0; i < num_rows; i++)
1157  {
1158  if (x_data[i] != 0)
1159  {
1160  for (jj = A_i[i]; jj < A_i[i+1]; jj++)
1161  {
1162  j = A_j[jj];
1163  /* y_data[j] += A_data[jj] * x_data[i]; */
1164  y_data[j] = 1;
1165  }
1166  }
1167  }
1168 }
1169 
1170 /* Based on hypre_ParCSRCommHandleCreate in hypre's par_csr_communication.c. The
1171  input variable job controls the communication type: 1=Matvec, 2=MatvecT. */
1172 hypre_ParCSRCommHandle *
1173 hypre_ParCSRCommHandleCreate_bool(HYPRE_Int job,
1174  hypre_ParCSRCommPkg *comm_pkg,
1175  HYPRE_Bool *send_data,
1176  HYPRE_Bool *recv_data)
1177 {
1178  HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1179  HYPRE_Int num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg);
1180  MPI_Comm comm = hypre_ParCSRCommPkgComm(comm_pkg);
1181 
1182  hypre_ParCSRCommHandle *comm_handle;
1183  HYPRE_Int num_requests;
1184  hypre_MPI_Request *requests;
1185 
1186  HYPRE_Int i, j;
1187  HYPRE_Int my_id, num_procs;
1188  HYPRE_Int ip, vec_start, vec_len;
1189 
1190  num_requests = num_sends + num_recvs;
1191  requests = mfem_hypre_CTAlloc(hypre_MPI_Request, num_requests);
1192 
1193  hypre_MPI_Comm_size(comm, &num_procs);
1194  hypre_MPI_Comm_rank(comm, &my_id);
1195 
1196  j = 0;
1197  switch (job)
1198  {
1199  case 1:
1200  {
1201  HYPRE_Bool *d_send_data = (HYPRE_Bool *) send_data;
1202  HYPRE_Bool *d_recv_data = (HYPRE_Bool *) recv_data;
1203  for (i = 0; i < num_recvs; i++)
1204  {
1205  ip = hypre_ParCSRCommPkgRecvProc(comm_pkg, i);
1206  vec_start = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i);
1207  vec_len = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i+1)-vec_start;
1208  hypre_MPI_Irecv(&d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1209  ip, 0, comm, &requests[j++]);
1210  }
1211  for (i = 0; i < num_sends; i++)
1212  {
1213  vec_start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1214  vec_len = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1)-vec_start;
1215  ip = hypre_ParCSRCommPkgSendProc(comm_pkg, i);
1216  hypre_MPI_Isend(&d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1217  ip, 0, comm, &requests[j++]);
1218  }
1219  break;
1220  }
1221  case 2:
1222  {
1223  HYPRE_Bool *d_send_data = (HYPRE_Bool *) send_data;
1224  HYPRE_Bool *d_recv_data = (HYPRE_Bool *) recv_data;
1225  for (i = 0; i < num_sends; i++)
1226  {
1227  vec_start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1228  vec_len = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1)-vec_start;
1229  ip = hypre_ParCSRCommPkgSendProc(comm_pkg, i);
1230  hypre_MPI_Irecv(&d_recv_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1231  ip, 0, comm, &requests[j++]);
1232  }
1233  for (i = 0; i < num_recvs; i++)
1234  {
1235  ip = hypre_ParCSRCommPkgRecvProc(comm_pkg, i);
1236  vec_start = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i);
1237  vec_len = hypre_ParCSRCommPkgRecvVecStart(comm_pkg, i+1)-vec_start;
1238  hypre_MPI_Isend(&d_send_data[vec_start], vec_len, HYPRE_MPI_BOOL,
1239  ip, 0, comm, &requests[j++]);
1240  }
1241  break;
1242  }
1243  }
1244  /*--------------------------------------------------------------------
1245  * set up comm_handle and return
1246  *--------------------------------------------------------------------*/
1247 
1248  comm_handle = mfem_hypre_CTAlloc(hypre_ParCSRCommHandle, 1);
1249 
1250  hypre_ParCSRCommHandleCommPkg(comm_handle) = comm_pkg;
1251  hypre_ParCSRCommHandleSendData(comm_handle) = send_data;
1252  hypre_ParCSRCommHandleRecvData(comm_handle) = recv_data;
1253  hypre_ParCSRCommHandleNumRequests(comm_handle) = num_requests;
1254  hypre_ParCSRCommHandleRequests(comm_handle) = requests;
1255 
1256  return comm_handle;
1257 }
1258 
1259 /* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
1260 void hypre_ParCSRMatrixBooleanMatvec(hypre_ParCSRMatrix *A,
1261  HYPRE_Bool alpha,
1262  HYPRE_Bool *x,
1263  HYPRE_Bool beta,
1264  HYPRE_Bool *y)
1265 {
1266  hypre_ParCSRCommHandle *comm_handle;
1267  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1268  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1269  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1270 
1271  HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1272  HYPRE_Int num_sends, i, j, index;
1273 
1274  HYPRE_Bool *x_tmp, *x_buf;
1275 
1276  x_tmp = mfem_hypre_CTAlloc(HYPRE_Bool, num_cols_offd);
1277 
1278  /*---------------------------------------------------------------------
1279  * If there exists no CommPkg for A, a CommPkg is generated using
1280  * equally load balanced partitionings
1281  *--------------------------------------------------------------------*/
1282  if (!comm_pkg)
1283  {
1284  hypre_MatvecCommPkgCreate(A);
1285  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1286  }
1287 
1288  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1289  x_buf = mfem_hypre_CTAlloc(
1290  HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1291 
1292  index = 0;
1293  for (i = 0; i < num_sends; i++)
1294  {
1295  j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
1296  for ( ; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
1297  {
1298  x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j)];
1299  }
1300  }
1301 
1302  comm_handle = hypre_ParCSRCommHandleCreate_bool(1, comm_pkg, x_buf, x_tmp);
1303 
1304  hypre_CSRMatrixBooleanMatvec(diag, alpha, x, beta, y);
1305 
1306  hypre_ParCSRCommHandleDestroy(comm_handle);
1307 
1308  if (num_cols_offd)
1309  {
1310  hypre_CSRMatrixBooleanMatvec(offd, alpha, x_tmp, 1, y);
1311  }
1312 
1313  mfem_hypre_TFree(x_buf);
1314  mfem_hypre_TFree(x_tmp);
1315 }
1316 
1317 /* Based on hypre_ParCSRMatrixMatvecT in par_csr_matvec.c */
1318 void hypre_ParCSRMatrixBooleanMatvecT(hypre_ParCSRMatrix *A,
1319  HYPRE_Bool alpha,
1320  HYPRE_Bool *x,
1321  HYPRE_Bool beta,
1322  HYPRE_Bool *y)
1323 {
1324  hypre_ParCSRCommHandle *comm_handle;
1325  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1326  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
1327  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
1328  HYPRE_Bool *y_tmp;
1329  HYPRE_Bool *y_buf;
1330 
1331  HYPRE_Int num_cols_offd = hypre_CSRMatrixNumCols(offd);
1332 
1333  HYPRE_Int i, j, jj, end, num_sends;
1334 
1335  y_tmp = mfem_hypre_TAlloc(HYPRE_Bool, num_cols_offd);
1336 
1337  /*---------------------------------------------------------------------
1338  * If there exists no CommPkg for A, a CommPkg is generated using
1339  * equally load balanced partitionings
1340  *--------------------------------------------------------------------*/
1341  if (!comm_pkg)
1342  {
1343  hypre_MatvecCommPkgCreate(A);
1344  comm_pkg = hypre_ParCSRMatrixCommPkg(A);
1345  }
1346 
1347  num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
1348  y_buf = mfem_hypre_CTAlloc(
1349  HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
1350 
1351  if (num_cols_offd)
1352  {
1353 #if MFEM_HYPRE_VERSION >= 21100
1354  if (A->offdT)
1355  {
1356  // offdT is optional. Used only if it's present.
1357  hypre_CSRMatrixBooleanMatvec(A->offdT, alpha, x, 0, y_tmp);
1358  }
1359  else
1360 #endif
1361  {
1362  hypre_CSRMatrixBooleanMatvecT(offd, alpha, x, 0, y_tmp);
1363  }
1364  }
1365 
1366  comm_handle = hypre_ParCSRCommHandleCreate_bool(2, comm_pkg, y_tmp, y_buf);
1367 
1368 #if MFEM_HYPRE_VERSION >= 21100
1369  if (A->diagT)
1370  {
1371  // diagT is optional. Used only if it's present.
1372  hypre_CSRMatrixBooleanMatvec(A->diagT, alpha, x, beta, y);
1373  }
1374  else
1375 #endif
1376  {
1377  hypre_CSRMatrixBooleanMatvecT(diag, alpha, x, beta, y);
1378  }
1379 
1380  hypre_ParCSRCommHandleDestroy(comm_handle);
1381 
1382  for (i = 0; i < num_sends; i++)
1383  {
1384  end = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1);
1385  for (j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i); j < end; j++)
1386  {
1387  jj = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
1388  y[jj] = y[jj] || y_buf[j];
1389  }
1390  }
1391 
1392  mfem_hypre_TFree(y_buf);
1393  mfem_hypre_TFree(y_tmp);
1394 }
1395 
1396 HYPRE_Int
1397 hypre_CSRMatrixSum(hypre_CSRMatrix *A,
1398  HYPRE_Complex beta,
1399  hypre_CSRMatrix *B)
1400 {
1401  HYPRE_Complex *A_data = hypre_CSRMatrixData(A);
1402  HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1403  HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
1404  HYPRE_Int nrows_A = hypre_CSRMatrixNumRows(A);
1405  HYPRE_Int ncols_A = hypre_CSRMatrixNumCols(A);
1406  HYPRE_Complex *B_data = hypre_CSRMatrixData(B);
1407  HYPRE_Int *B_i = hypre_CSRMatrixI(B);
1408  HYPRE_Int *B_j = hypre_CSRMatrixJ(B);
1409  HYPRE_Int nrows_B = hypre_CSRMatrixNumRows(B);
1410  HYPRE_Int ncols_B = hypre_CSRMatrixNumCols(B);
1411 
1412  HYPRE_Int ia, j, pos;
1413  HYPRE_Int *marker;
1414 
1415  if (nrows_A != nrows_B || ncols_A != ncols_B)
1416  {
1417  return -1; /* error: incompatible matrix dimensions */
1418  }
1419 
1420  marker = mfem_hypre_CTAlloc(HYPRE_Int, ncols_A);
1421  for (ia = 0; ia < ncols_A; ia++)
1422  {
1423  marker[ia] = -1;
1424  }
1425 
1426  for (ia = 0; ia < nrows_A; ia++)
1427  {
1428  for (j = A_i[ia]; j < A_i[ia+1]; j++)
1429  {
1430  marker[A_j[j]] = j;
1431  }
1432 
1433  for (j = B_i[ia]; j < B_i[ia+1]; j++)
1434  {
1435  pos = marker[B_j[j]];
1436  if (pos < A_i[ia])
1437  {
1438  return -2; /* error: found an entry in B that is not present in A */
1439  }
1440  A_data[pos] += beta * B_data[j];
1441  }
1442  }
1443 
1444  mfem_hypre_TFree(marker);
1445  return 0;
1446 }
1447 
1448 hypre_ParCSRMatrix *
1449 hypre_ParCSRMatrixAdd(hypre_ParCSRMatrix *A,
1450  hypre_ParCSRMatrix *B)
1451 {
1452  MPI_Comm comm = hypre_ParCSRMatrixComm(A);
1453  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1454  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
1455  HYPRE_Int *A_cmap = hypre_ParCSRMatrixColMapOffd(A);
1456  HYPRE_Int A_cmap_size = hypre_CSRMatrixNumCols(A_offd);
1457  hypre_CSRMatrix *B_diag = hypre_ParCSRMatrixDiag(B);
1458  hypre_CSRMatrix *B_offd = hypre_ParCSRMatrixOffd(B);
1459  HYPRE_Int *B_cmap = hypre_ParCSRMatrixColMapOffd(B);
1460  HYPRE_Int B_cmap_size = hypre_CSRMatrixNumCols(B_offd);
1461  hypre_ParCSRMatrix *C;
1462  hypre_CSRMatrix *C_diag;
1463  hypre_CSRMatrix *C_offd;
1464  HYPRE_Int *C_cmap;
1465  HYPRE_Int im;
1466  HYPRE_Int cmap_differ;
1467 
1468  /* Check if A_cmap and B_cmap are the same. */
1469  cmap_differ = 0;
1470  if (A_cmap_size != B_cmap_size)
1471  {
1472  cmap_differ = 1; /* A and B have different cmap_size */
1473  }
1474  else
1475  {
1476  for (im = 0; im < A_cmap_size; im++)
1477  {
1478  if (A_cmap[im] != B_cmap[im])
1479  {
1480  cmap_differ = 1; /* A and B have different cmap arrays */
1481  break;
1482  }
1483  }
1484  }
1485 
1486  if ( cmap_differ == 0 )
1487  {
1488  /* A and B have the same column mapping for their off-diagonal blocks so
1489  we can sum the diagonal and off-diagonal blocks separately and reduce
1490  temporary memory usage. */
1491 
1492  /* Add diagonals, off-diagonals, copy cmap. */
1493  C_diag = hypre_CSRMatrixAdd(A_diag, B_diag);
1494  if (!C_diag)
1495  {
1496  return NULL; /* error: A_diag and B_diag have different dimensions */
1497  }
1498  C_offd = hypre_CSRMatrixAdd(A_offd, B_offd);
1499  if (!C_offd)
1500  {
1501  hypre_CSRMatrixDestroy(C_diag);
1502  return NULL; /* error: A_offd and B_offd have different dimensions */
1503  }
1504  /* copy A_cmap -> C_cmap */
1505  C_cmap = mfem_hypre_TAlloc(HYPRE_Int, A_cmap_size);
1506  for (im = 0; im < A_cmap_size; im++)
1507  {
1508  C_cmap[im] = A_cmap[im];
1509  }
1510 
1511  C = hypre_ParCSRMatrixCreate(comm,
1512  hypre_ParCSRMatrixGlobalNumRows(A),
1513  hypre_ParCSRMatrixGlobalNumCols(A),
1514  hypre_ParCSRMatrixRowStarts(A),
1515  hypre_ParCSRMatrixColStarts(A),
1516  hypre_CSRMatrixNumCols(C_offd),
1517  hypre_CSRMatrixNumNonzeros(C_diag),
1518  hypre_CSRMatrixNumNonzeros(C_offd));
1519 
1520  /* In C, destroy diag/offd (allocated by Create) and replace them with
1521  C_diag/C_offd. */
1522  hypre_CSRMatrixDestroy(hypre_ParCSRMatrixDiag(C));
1523  hypre_CSRMatrixDestroy(hypre_ParCSRMatrixOffd(C));
1524  hypre_ParCSRMatrixDiag(C) = C_diag;
1525  hypre_ParCSRMatrixOffd(C) = C_offd;
1526 
1527  hypre_ParCSRMatrixColMapOffd(C) = C_cmap;
1528  }
1529  else
1530  {
1531  /* A and B have different column mappings for their off-diagonal blocks so
1532  we need to use the column maps to create full-width CSR matrices. */
1533 
1534  int ierr = 0;
1535  hypre_CSRMatrix * csr_A;
1536  hypre_CSRMatrix * csr_B;
1537  hypre_CSRMatrix * csr_C_temp;
1538 
1539  /* merge diag and off-diag portions of A */
1540  csr_A = hypre_MergeDiagAndOffd(A);
1541 
1542  /* merge diag and off-diag portions of B */
1543  csr_B = hypre_MergeDiagAndOffd(B);
1544 
1545  /* add A and B */
1546  csr_C_temp = hypre_CSRMatrixAdd(csr_A,csr_B);
1547 
1548  /* delete CSR versions of A and B */
1549  ierr += hypre_CSRMatrixDestroy(csr_A);
1550  ierr += hypre_CSRMatrixDestroy(csr_B);
1551 
1552  /* create a new empty ParCSR matrix to contain the sum */
1553  C = hypre_ParCSRMatrixCreate(hypre_ParCSRMatrixComm(A),
1554  hypre_ParCSRMatrixGlobalNumRows(A),
1555  hypre_ParCSRMatrixGlobalNumCols(A),
1556  hypre_ParCSRMatrixRowStarts(A),
1557  hypre_ParCSRMatrixColStarts(A),
1558  0, 0, 0);
1559 
1560  /* split C into diag and off-diag portions */
1561  /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the
1562  number of columns in csr_C_temp which is the global number of columns
1563  in A and B. This does not scale well. */
1564  ierr += GenerateDiagAndOffd(csr_C_temp, C,
1565  hypre_ParCSRMatrixFirstColDiag(A),
1566  hypre_ParCSRMatrixLastColDiag(A));
1567 
1568  /* delete CSR version of C */
1569  ierr += hypre_CSRMatrixDestroy(csr_C_temp);
1570 
1571  MFEM_VERIFY(ierr == 0, "");
1572  }
1573 
1574  /* hypre_ParCSRMatrixSetNumNonzeros(A); */
1575 
1576  /* Make sure that the first entry in each row is the diagonal one. */
1577  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(C));
1578 
1579  /* C owns diag, offd, and cmap. */
1580  hypre_ParCSRMatrixSetDataOwner(C, 1);
1581  /* C does not own row and column starts. */
1582  hypre_ParCSRMatrixSetRowStartsOwner(C, 0);
1583  hypre_ParCSRMatrixSetColStartsOwner(C, 0);
1584 
1585  return C;
1586 }
1587 
1588 HYPRE_Int
1589 hypre_ParCSRMatrixSum(hypre_ParCSRMatrix *A,
1590  HYPRE_Complex beta,
1591  hypre_ParCSRMatrix *B)
1592 {
1593  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1594  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
1595  hypre_CSRMatrix *B_diag = hypre_ParCSRMatrixDiag(B);
1596  hypre_CSRMatrix *B_offd = hypre_ParCSRMatrixOffd(B);
1597  HYPRE_Int ncols_B_offd = hypre_CSRMatrixNumCols(B_offd);
1598  HYPRE_Int error;
1599 
1600  error = hypre_CSRMatrixSum(A_diag, beta, B_diag);
1601  if (ncols_B_offd > 0) /* treat B_offd as zero if it has no columns */
1602  {
1603  error = error ? error : hypre_CSRMatrixSum(A_offd, beta, B_offd);
1604  }
1605 
1606  return error;
1607 }
1608 
1609 HYPRE_Int
1610 hypre_CSRMatrixSetConstantValues(hypre_CSRMatrix *A,
1611  HYPRE_Complex value)
1612 {
1613  HYPRE_Complex *A_data = hypre_CSRMatrixData(A);
1614  HYPRE_Int A_nnz = hypre_CSRMatrixNumNonzeros(A);
1615  HYPRE_Int ia;
1616 
1617  for (ia = 0; ia < A_nnz; ia++)
1618  {
1619  A_data[ia] = value;
1620  }
1621 
1622  return 0;
1623 }
1624 
1625 HYPRE_Int
1626 hypre_ParCSRMatrixSetConstantValues(hypre_ParCSRMatrix *A,
1627  HYPRE_Complex value)
1628 {
1629  hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixDiag(A), value);
1630  hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixOffd(A), value);
1631 
1632  return 0;
1633 }
1634 
1635 } // namespace mfem::internal
1636 
1637 } // namespace mfem
1638 
1639 #endif // MFEM_USE_MPI
hypre_parcsr.hpp
a
double a
Definition: lissajous.cpp:41
index
int index(int i, int j, int nx, int ny)
Definition: life.cpp:241
alpha
const double alpha
Definition: ex15.cpp:336
Generated on Fri Feb 18 2022 18:54:18 for MFEM by   doxygen 1.8.5