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