MFEM  v4.6.0
Finite element discretization library
complex_operator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 "complex_operator.hpp"
13 #include <set>
14 #include <map>
15 
16 namespace mfem
17 {
18 
20  bool ownReal, bool ownImag,
21  Convention convention)
22  : Operator(2*((Op_Real)?Op_Real->Height():Op_Imag->Height()),
23  2*((Op_Real)?Op_Real->Width():Op_Imag->Width()))
24  , Op_Real_(Op_Real)
25  , Op_Imag_(Op_Imag)
26  , ownReal_(ownReal)
27  , ownImag_(ownImag)
28  , convention_(convention)
29  , x_r_()
30  , x_i_()
31  , y_r_()
32  , y_i_()
33  , u_(NULL)
34  , v_(NULL)
35 {}
36 
38 {
39  if (ownReal_) { delete Op_Real_; }
40  if (ownImag_) { delete Op_Imag_; }
41  delete u_;
42  delete v_;
43 }
44 
46 {
47  MFEM_ASSERT(Op_Real_, "ComplexOperator has no real part!");
48  return *Op_Real_;
49 }
50 
52 {
53  MFEM_ASSERT(Op_Imag_, "ComplexOperator has no imaginary part!");
54  return *Op_Imag_;
55 }
56 
58 {
59  MFEM_ASSERT(Op_Real_, "ComplexOperator has no real part!");
60  return *Op_Real_;
61 }
62 
64 {
65  MFEM_ASSERT(Op_Imag_, "ComplexOperator has no imaginary part!");
66  return *Op_Imag_;
67 }
68 
69 void ComplexOperator::Mult(const Vector &x, Vector &y) const
70 {
71  x.Read();
72  y.UseDevice(true); y = 0.0;
73 
74  x_r_.MakeRef(const_cast<Vector&>(x), 0, width/2);
75  x_i_.MakeRef(const_cast<Vector&>(x), width/2, width/2);
76 
77  y_r_.MakeRef(y, 0, height/2);
78  y_i_.MakeRef(y, height/2, height/2);
79 
80  this->Mult(x_r_, x_i_, y_r_, y_i_);
81 
84 }
85 
86 void ComplexOperator::Mult(const Vector &x_r, const Vector &x_i,
87  Vector &y_r, Vector &y_i) const
88 {
89  if (Op_Real_)
90  {
91  Op_Real_->Mult(x_r, y_r);
92  Op_Real_->Mult(x_i, y_i);
93  }
94  else
95  {
96  y_r = 0.0;
97  y_i = 0.0;
98  }
99 
100  if (Op_Imag_)
101  {
102  if (!v_) { v_ = new Vector(); }
103  v_->UseDevice(true);
104  v_->SetSize(Op_Imag_->Height());
105 
106  Op_Imag_->Mult(x_i, *v_);
107  y_r.Add(-1.0, *v_);
108  Op_Imag_->Mult(x_r, *v_);
109  y_i.Add(1.0, *v_);
110  }
111 
113  {
114  y_i *= -1.0;
115  }
116 }
117 
119 {
120  x.Read();
121  y.UseDevice(true); y = 0.0;
122 
123  x_r_.MakeRef(const_cast<Vector&>(x), 0, height/2);
124  x_i_.MakeRef(const_cast<Vector&>(x), height/2, height/2);
125 
126  y_r_.MakeRef(y, 0, width/2);
127  y_i_.MakeRef(y, width/2, width/2);
128 
129  this->MultTranspose(x_r_, x_i_, y_r_, y_i_);
130 
133 }
134 
135 void ComplexOperator::MultTranspose(const Vector &x_r, const Vector &x_i,
136  Vector &y_r, Vector &y_i) const
137 {
138  if (Op_Real_)
139  {
140  Op_Real_->MultTranspose(x_r, y_r);
141  Op_Real_->MultTranspose(x_i, y_i);
142 
144  {
145  y_i *= -1.0;
146  }
147  }
148  else
149  {
150  y_r = 0.0;
151  y_i = 0.0;
152  }
153 
154  if (Op_Imag_)
155  {
156  if (!u_) { u_ = new Vector(); }
157  u_->UseDevice(true);
158  u_->SetSize(Op_Imag_->Width());
159 
160  Op_Imag_->MultTranspose(x_i, *u_);
161  y_r.Add(convention_ == BLOCK_SYMMETRIC ? -1.0 : 1.0, *u_);
162  Op_Imag_->MultTranspose(x_r, *u_);
163  y_i.Add(-1.0, *u_);
164  }
165 }
166 
167 
169 {
170  MFEM_ASSERT(Op_Real_, "ComplexSparseMatrix has no real part!");
171  return dynamic_cast<SparseMatrix &>(*Op_Real_);
172 }
173 
175 {
176  MFEM_ASSERT(Op_Imag_, "ComplexSparseMatrix has no imaginary part!");
177  return dynamic_cast<SparseMatrix &>(*Op_Imag_);
178 }
179 
181 {
182  MFEM_ASSERT(Op_Real_, "ComplexSparseMatrix has no real part!");
183  return dynamic_cast<const SparseMatrix &>(*Op_Real_);
184 }
185 
187 {
188  MFEM_ASSERT(Op_Imag_, "ComplexSparseMatrix has no imaginary part!");
189  return dynamic_cast<const SparseMatrix &>(*Op_Imag_);
190 }
191 
193 {
194  SparseMatrix * A_r = dynamic_cast<SparseMatrix*>(Op_Real_);
195  SparseMatrix * A_i = dynamic_cast<SparseMatrix*>(Op_Imag_);
196 
197  const int nrows_r = (A_r)?A_r->Height():0;
198  const int nrows_i = (A_i)?A_i->Height():0;
199  const int nrows = std::max(nrows_r, nrows_i);
200 
201  const int *I_r = (A_r)?A_r->GetI():NULL;
202  const int *I_i = (A_i)?A_i->GetI():NULL;
203 
204  const int *J_r = (A_r)?A_r->GetJ():NULL;
205  const int *J_i = (A_i)?A_i->GetJ():NULL;
206 
207  const double *D_r = (A_r)?A_r->GetData():NULL;
208  const double *D_i = (A_i)?A_i->GetData():NULL;
209 
210  const int nnz_r = (I_r)?I_r[nrows]:0;
211  const int nnz_i = (I_i)?I_i[nrows]:0;
212  const int nnz = 2 * (nnz_r + nnz_i);
213 
214  int *I = Memory<int>(this->Height()+1);
215  int *J = Memory<int>(nnz);
216  double *D = Memory<double>(nnz);
217 
218  const double factor = (convention_ == HERMITIAN) ? 1.0 : -1.0;
219 
220  I[0] = 0;
221  I[nrows] = nnz_r + nnz_i;
222  for (int i=0; i<nrows; i++)
223  {
224  I[i + 1] = ((I_r)?I_r[i+1]:0) + ((I_i)?I_i[i+1]:0);
225  I[i + nrows + 1] = I[i+1] + nnz_r + nnz_i;
226 
227  if (I_r)
228  {
229  const int off_i = (I_i)?(I_i[i+1] - I_i[i]):0;
230  for (int j=0; j<I_r[i+1] - I_r[i]; j++)
231  {
232  J[I[i] + j] = J_r[I_r[i] + j];
233  D[I[i] + j] = D_r[I_r[i] + j];
234 
235  J[I[i+nrows] + off_i + j] = J_r[I_r[i] + j] + nrows;
236  D[I[i+nrows] + off_i + j] = factor*D_r[I_r[i] + j];
237  }
238  }
239  if (I_i)
240  {
241  const int off_r = (I_r)?(I_r[i+1] - I_r[i]):0;
242  for (int j=0; j<I_i[i+1] - I_i[i]; j++)
243  {
244  J[I[i] + off_r + j] = J_i[I_i[i] + j] + nrows;
245  D[I[i] + off_r + j] = -D_i[I_i[i] + j];
246 
247  J[I[i+nrows] + j] = J_i[I_i[i] + j];
248  D[I[i+nrows] + j] = factor*D_i[I_i[i] + j];
249  }
250  }
251  }
252 
253  return new SparseMatrix(I, J, D, this->Height(), this->Width());
254 }
255 
256 
257 #ifdef MFEM_USE_SUITESPARSE
258 
260 {
261  mat = NULL;
262  Numeric = NULL;
263  AI = AJ = NULL;
264  if (!use_long_ints)
265  {
266  umfpack_zi_defaults(Control);
267  }
268  else
269  {
270  umfpack_zl_defaults(Control);
271  }
272 }
273 
275 {
276  void *Symbolic;
277 
278  if (Numeric)
279  {
280  if (!use_long_ints)
281  {
282  umfpack_zi_free_numeric(&Numeric);
283  }
284  else
285  {
286  umfpack_zl_free_numeric(&Numeric);
287  }
288  }
289 
290  mat = const_cast<ComplexSparseMatrix *>
291  (dynamic_cast<const ComplexSparseMatrix *>(&op));
292  MFEM_VERIFY(mat, "not a ComplexSparseMatrix");
293 
294  MFEM_VERIFY(mat->real().NumNonZeroElems() == mat->imag().NumNonZeroElems(),
295  "Real and imag Sparsity pattern mismatch: Try setting Assemble (skip_zeros = 0)");
296 
297  // UMFPack requires that the column-indices in mat corresponding to each
298  // row be sorted.
299  // Generally, this will modify the ordering of the entries of mat.
300 
303 
304  height = mat->real().Height();
305  width = mat->real().Width();
306  MFEM_VERIFY(width == height, "not a square matrix");
307 
308  const int * Ap =
309  mat->real().HostReadI(); // assuming real and imag have the same sparsity
310  const int * Ai = mat->real().HostReadJ();
311  const double * Ax = mat->real().HostReadData();
312  const double * Az = mat->imag().HostReadData();
313 
314  if (!use_long_ints)
315  {
316  int status = umfpack_zi_symbolic(width,width,Ap,Ai,Ax,Az,&Symbolic,
317  Control,Info);
318  if (status < 0)
319  {
320  umfpack_zi_report_info(Control, Info);
321  umfpack_zi_report_status(Control, status);
322  mfem_error("ComplexUMFPackSolver::SetOperator :"
323  " umfpack_zi_symbolic() failed!");
324  }
325 
326  status = umfpack_zi_numeric(Ap, Ai, Ax, Az, Symbolic, &Numeric,
327  Control, Info);
328  if (status < 0)
329  {
330  umfpack_zi_report_info(Control, Info);
331  umfpack_zi_report_status(Control, status);
332  mfem_error("ComplexUMFPackSolver::SetOperator :"
333  " umfpack_zi_numeric() failed!");
334  }
335  umfpack_zi_free_symbolic(&Symbolic);
336  }
337  else
338  {
339  SuiteSparse_long status;
340 
341  delete [] AJ;
342  delete [] AI;
343  AI = new SuiteSparse_long[width + 1];
344  AJ = new SuiteSparse_long[Ap[width]];
345  for (int i = 0; i <= width; i++)
346  {
347  AI[i] = (SuiteSparse_long)(Ap[i]);
348  }
349  for (int i = 0; i < Ap[width]; i++)
350  {
351  AJ[i] = (SuiteSparse_long)(Ai[i]);
352  }
353 
354  status = umfpack_zl_symbolic(width, width, AI, AJ, Ax, Az, &Symbolic,
355  Control, Info);
356  if (status < 0)
357  {
358  umfpack_zl_report_info(Control, Info);
359  umfpack_zl_report_status(Control, status);
360  mfem_error("ComplexUMFPackSolver::SetOperator :"
361  " umfpack_zl_symbolic() failed!");
362  }
363 
364  status = umfpack_zl_numeric(AI, AJ, Ax, Az, Symbolic, &Numeric,
365  Control, Info);
366  if (status < 0)
367  {
368  umfpack_zl_report_info(Control, Info);
369  umfpack_zl_report_status(Control, status);
370  mfem_error("ComplexUMFPackSolver::SetOperator :"
371  " umfpack_zl_numeric() failed!");
372  }
373  umfpack_zl_free_symbolic(&Symbolic);
374  }
375 }
376 
378 {
379  if (mat == NULL)
380  mfem_error("ComplexUMFPackSolver::Mult : matrix is not set!"
381  " Call SetOperator first!");
382 
383  b.HostRead();
384  x.HostReadWrite();
385 
386  int n = b.Size()/2;
387  double * datax = x.GetData();
388  double * datab = b.GetData();
389 
390  // For the Block Symmetric case data the imaginary part
391  // has to be scaled by -1
393  Vector bimag;
394  if (conv == ComplexOperator::Convention::BLOCK_SYMMETRIC)
395  {
396  bimag.SetDataAndSize(&datab[n],n);
397  bimag *=-1.0;
398  }
399 
400  // Solve the transpose, since UMFPack expects CCS instead of CRS format
401  if (!use_long_ints)
402  {
403  int status =
404  umfpack_zi_solve(UMFPACK_Aat, mat->real().HostReadI(), mat->real().HostReadJ(),
406  datax, &datax[n], datab, &datab[n], Numeric, Control, Info);
407  umfpack_zi_report_info(Control, Info);
408  if (status < 0)
409  {
410  umfpack_zi_report_status(Control, status);
411  mfem_error("ComplexUMFPackSolver::Mult : umfpack_zi_solve() failed!");
412  }
413  }
414  else
415  {
416  SuiteSparse_long status =
417  umfpack_zl_solve(UMFPACK_Aat,AI,AJ,mat->real().HostReadData(),
418  mat->imag().HostReadData(),
419  datax,&datax[n],datab,&datab[n],Numeric,Control,Info);
420 
421  umfpack_zl_report_info(Control, Info);
422  if (status < 0)
423  {
424  umfpack_zl_report_status(Control, status);
425  mfem_error("ComplexUMFPackSolver::Mult : umfpack_zl_solve() failed!");
426  }
427  }
428  if (conv == ComplexOperator::Convention::BLOCK_SYMMETRIC)
429  {
430  bimag *=-1.0;
431  }
432 }
433 
435 {
436  if (mat == NULL)
437  mfem_error("ComplexUMFPackSolver::Mult : matrix is not set!"
438  " Call SetOperator first!");
439  b.HostRead();
440  x.HostReadWrite();
441  int n = b.Size()/2;
442  double * datax = x.GetData();
443  double * datab = b.GetData();
444 
446  Vector bimag;
447  bimag.SetDataAndSize(&datab[n],n);
448 
449  // Solve the Adjoint A^H x = b by solving
450  // the conjugate problem A^T \bar{x} = \bar{b}
451  if ((!transa && conv == ComplexOperator::HERMITIAN) ||
453  {
454  bimag *=-1.0;
455  }
456 
457  if (!use_long_ints)
458  {
459  int status =
460  umfpack_zi_solve(UMFPACK_A, mat->real().HostReadI(), mat->real().HostReadJ(),
462  datax, &datax[n], datab, &datab[n], Numeric, Control, Info);
463  umfpack_zi_report_info(Control, Info);
464  if (status < 0)
465  {
466  umfpack_zi_report_status(Control, status);
467  mfem_error("ComplexUMFPackSolver::Mult : umfpack_zi_solve() failed!");
468  }
469  }
470  else
471  {
472  SuiteSparse_long status =
473  umfpack_zl_solve(UMFPACK_A,AI,AJ,mat->real().HostReadData(),
474  mat->imag().HostReadData(),
475  datax,&datax[n],datab,&datab[n],Numeric,Control,Info);
476 
477  umfpack_zl_report_info(Control, Info);
478  if (status < 0)
479  {
480  umfpack_zl_report_status(Control, status);
481  mfem_error("ComplexUMFPackSolver::Mult : umfpack_zl_solve() failed!");
482  }
483  }
484  if (!transa)
485  {
486  Vector ximag;
487  ximag.SetDataAndSize(&datax[n],n);
488  ximag *=-1.0;
489  }
490  if ((!transa && conv == ComplexOperator::HERMITIAN) ||
492  {
493  bimag *=-1.0;
494  }
495 }
496 
498 {
499  delete [] AJ;
500  delete [] AI;
501  if (Numeric)
502  {
503  if (!use_long_ints)
504  {
505  umfpack_zi_free_numeric(&Numeric);
506  }
507  else
508  {
509  umfpack_zl_free_numeric(&Numeric);
510  }
511  }
512 }
513 
514 #endif
515 
516 #ifdef MFEM_USE_MPI
517 
519  HypreParMatrix * A_Imag,
520  bool ownReal, bool ownImag,
521  Convention convention)
522  : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
523 {
524  comm_ = (A_Real) ? A_Real->GetComm() :
525  ((A_Imag) ? A_Imag->GetComm() : MPI_COMM_WORLD);
526 
527  MPI_Comm_rank(comm_, &myid_);
528  MPI_Comm_size(comm_, &nranks_);
529 }
530 
532 {
533  MFEM_ASSERT(Op_Real_, "ComplexHypreParMatrix has no real part!");
534  return dynamic_cast<HypreParMatrix &>(*Op_Real_);
535 }
536 
538 {
539  MFEM_ASSERT(Op_Imag_, "ComplexHypreParMatrix has no imaginary part!");
540  return dynamic_cast<HypreParMatrix &>(*Op_Imag_);
541 }
542 
544 {
545  MFEM_ASSERT(Op_Real_, "ComplexHypreParMatrix has no real part!");
546  return dynamic_cast<const HypreParMatrix &>(*Op_Real_);
547 }
548 
550 {
551  MFEM_ASSERT(Op_Imag_, "ComplexHypreParMatrix has no imaginary part!");
552  return dynamic_cast<const HypreParMatrix &>(*Op_Imag_);
553 }
554 
556 {
557  HypreParMatrix * A_r = dynamic_cast<HypreParMatrix*>(Op_Real_);
558  HypreParMatrix * A_i = dynamic_cast<HypreParMatrix*>(Op_Imag_);
559 
560  if ( A_r == NULL && A_i == NULL ) { return NULL; }
561 
562  HYPRE_BigInt global_num_rows_r = (A_r) ? A_r->GetGlobalNumRows() : 0;
563  HYPRE_BigInt global_num_rows_i = (A_i) ? A_i->GetGlobalNumRows() : 0;
564  HYPRE_BigInt global_num_rows = std::max(global_num_rows_r,
565  global_num_rows_i);
566 
567  HYPRE_BigInt global_num_cols_r = (A_r) ? A_r->GetGlobalNumCols() : 0;
568  HYPRE_BigInt global_num_cols_i = (A_i) ? A_i->GetGlobalNumCols() : 0;
569  HYPRE_BigInt global_num_cols = std::max(global_num_cols_r,
570  global_num_cols_i);
571 
572  int row_starts_size = (HYPRE_AssumedPartitionCheck()) ? 2 : nranks_ + 1;
573  HYPRE_BigInt * row_starts = mfem_hypre_CTAlloc_host(HYPRE_BigInt,
574  row_starts_size);
575  HYPRE_BigInt * col_starts = mfem_hypre_CTAlloc_host(HYPRE_BigInt,
576  row_starts_size);
577 
578  const HYPRE_BigInt * row_starts_z = (A_r) ? A_r->RowPart() :
579  ((A_i) ? A_i->RowPart() : NULL);
580  const HYPRE_BigInt * col_starts_z = (A_r) ? A_r->ColPart() :
581  ((A_i) ? A_i->ColPart() : NULL);
582 
583  for (int i = 0; i < row_starts_size; i++)
584  {
585  row_starts[i] = 2 * row_starts_z[i];
586  col_starts[i] = 2 * col_starts_z[i];
587  }
588 
589  SparseMatrix diag_r, diag_i, offd_r, offd_i;
590  HYPRE_BigInt * cmap_r = NULL, * cmap_i = NULL;
591 
592  int nrows_r = 0, nrows_i = 0, ncols_r = 0, ncols_i = 0;
593  int ncols_offd_r = 0, ncols_offd_i = 0;
594  if (A_r)
595  {
596  A_r->GetDiag(diag_r);
597  A_r->GetOffd(offd_r, cmap_r);
598  nrows_r = diag_r.Height();
599  ncols_r = diag_r.Width();
600  ncols_offd_r = offd_r.Width();
601  }
602  if (A_i)
603  {
604  A_i->GetDiag(diag_i);
605  A_i->GetOffd(offd_i, cmap_i);
606  nrows_i = diag_i.Height();
607  ncols_i = diag_i.Width();
608  ncols_offd_i = offd_i.Width();
609  }
610  int nrows = std::max(nrows_r, nrows_i);
611  int ncols = std::max(ncols_r, ncols_i);
612 
613  // Determine the unique set of off-diagonal columns global indices
614  std::set<HYPRE_BigInt> cset;
615  for (int i=0; i<ncols_offd_r; i++)
616  {
617  cset.insert(cmap_r[i]);
618  }
619  for (int i=0; i<ncols_offd_i; i++)
620  {
621  cset.insert(cmap_i[i]);
622  }
623  int num_cols_offd = (int)cset.size();
624 
625  // Extract pointers to the various CSR arrays of the diagonal blocks
626  const int * diag_r_I = (A_r) ? diag_r.GetI() : NULL;
627  const int * diag_i_I = (A_i) ? diag_i.GetI() : NULL;
628 
629  const int * diag_r_J = (A_r) ? diag_r.GetJ() : NULL;
630  const int * diag_i_J = (A_i) ? diag_i.GetJ() : NULL;
631 
632  const double * diag_r_D = (A_r) ? diag_r.GetData() : NULL;
633  const double * diag_i_D = (A_i) ? diag_i.GetData() : NULL;
634 
635  int diag_r_nnz = (diag_r_I) ? diag_r_I[nrows] : 0;
636  int diag_i_nnz = (diag_i_I) ? diag_i_I[nrows] : 0;
637  int diag_nnz = 2 * (diag_r_nnz + diag_i_nnz);
638 
639  // Extract pointers to the various CSR arrays of the off-diagonal blocks
640  const int * offd_r_I = (A_r) ? offd_r.GetI() : NULL;
641  const int * offd_i_I = (A_i) ? offd_i.GetI() : NULL;
642 
643  const int * offd_r_J = (A_r) ? offd_r.GetJ() : NULL;
644  const int * offd_i_J = (A_i) ? offd_i.GetJ() : NULL;
645 
646  const double * offd_r_D = (A_r) ? offd_r.GetData() : NULL;
647  const double * offd_i_D = (A_i) ? offd_i.GetData() : NULL;
648 
649  int offd_r_nnz = (offd_r_I) ? offd_r_I[nrows] : 0;
650  int offd_i_nnz = (offd_i_I) ? offd_i_I[nrows] : 0;
651  int offd_nnz = 2 * (offd_r_nnz + offd_i_nnz);
652 
653  // Allocate CSR arrays for the combined matrix
654  HYPRE_Int * diag_I = mfem_hypre_CTAlloc_host(HYPRE_Int, 2 * nrows + 1);
655  HYPRE_Int * diag_J = mfem_hypre_CTAlloc_host(HYPRE_Int, diag_nnz);
656  double * diag_D = mfem_hypre_CTAlloc_host(double, diag_nnz);
657 
658  HYPRE_Int * offd_I = mfem_hypre_CTAlloc_host(HYPRE_Int, 2 * nrows + 1);
659  HYPRE_Int * offd_J = mfem_hypre_CTAlloc_host(HYPRE_Int, offd_nnz);
660  double * offd_D = mfem_hypre_CTAlloc_host(double, offd_nnz);
661  HYPRE_BigInt * cmap = mfem_hypre_CTAlloc_host(HYPRE_BigInt,
662  2 * num_cols_offd);
663 
664  // Fill the CSR arrays for the diagonal portion of the matrix
665  const double factor = (convention_ == HERMITIAN) ? 1.0 : -1.0;
666 
667  diag_I[0] = 0;
668  diag_I[nrows] = diag_r_nnz + diag_i_nnz;
669  for (int i=0; i<nrows; i++)
670  {
671  diag_I[i + 1] = ((diag_r_I)?diag_r_I[i+1]:0) +
672  ((diag_i_I)?diag_i_I[i+1]:0);
673  diag_I[i + nrows + 1] = diag_I[i+1] + diag_r_nnz + diag_i_nnz;
674 
675  if (diag_r_I)
676  {
677  for (int j=0; j<diag_r_I[i+1] - diag_r_I[i]; j++)
678  {
679  diag_J[diag_I[i] + j] = diag_r_J[diag_r_I[i] + j];
680  diag_D[diag_I[i] + j] = diag_r_D[diag_r_I[i] + j];
681 
682  diag_J[diag_I[i+nrows] + j] =
683  diag_r_J[diag_r_I[i] + j] + ncols;
684  diag_D[diag_I[i+nrows] + j] =
685  factor * diag_r_D[diag_r_I[i] + j];
686  }
687  }
688  if (diag_i_I)
689  {
690  const int off_r = (diag_r_I)?(diag_r_I[i+1] - diag_r_I[i]):0;
691  for (int j=0; j<diag_i_I[i+1] - diag_i_I[i]; j++)
692  {
693  diag_J[diag_I[i] + off_r + j] = diag_i_J[diag_i_I[i] + j] + ncols;
694  diag_D[diag_I[i] + off_r + j] = -diag_i_D[diag_i_I[i] + j];
695 
696  diag_J[diag_I[i+nrows] + off_r + j] = diag_i_J[diag_i_I[i] + j];
697  diag_D[diag_I[i+nrows] + off_r + j] =
698  factor * diag_i_D[diag_i_I[i] + j];
699  }
700  }
701  }
702 
703  // Determine the mappings describing the layout of off-diagonal columns
704  int num_recv_procs = 0;
705  HYPRE_BigInt * offd_col_start_stop = NULL;
706  this->getColStartStop(A_r, A_i, num_recv_procs, offd_col_start_stop);
707 
708  std::set<HYPRE_BigInt>::iterator sit;
709  std::map<HYPRE_BigInt,HYPRE_BigInt> cmapa, cmapb, cinvmap;
710  for (sit=cset.begin(); sit!=cset.end(); sit++)
711  {
712  HYPRE_BigInt col_orig = *sit;
713  HYPRE_BigInt col_2x2 = -1;
714  HYPRE_BigInt col_size = 0;
715  for (int i=0; i<num_recv_procs; i++)
716  {
717  if (offd_col_start_stop[2*i] <= col_orig &&
718  col_orig < offd_col_start_stop[2*i+1])
719  {
720  col_2x2 = offd_col_start_stop[2*i] + col_orig;
721  col_size = offd_col_start_stop[2*i+1] - offd_col_start_stop[2*i];
722  break;
723  }
724  }
725  cmapa[*sit] = col_2x2;
726  cmapb[*sit] = col_2x2 + col_size;
727  cinvmap[col_2x2] = -1;
728  cinvmap[col_2x2 + col_size] = -1;
729  }
730  delete [] offd_col_start_stop;
731 
732  {
733  std::map<HYPRE_BigInt, HYPRE_BigInt>::iterator mit;
734  HYPRE_BigInt i = 0;
735  for (mit=cinvmap.begin(); mit!=cinvmap.end(); mit++, i++)
736  {
737  mit->second = i;
738  cmap[i] = mit->first;
739  }
740  }
741 
742  // Fill the CSR arrays for the off-diagonal portion of the matrix
743  offd_I[0] = 0;
744  offd_I[nrows] = offd_r_nnz + offd_i_nnz;
745  for (int i=0; i<nrows; i++)
746  {
747  offd_I[i + 1] = ((offd_r_I)?offd_r_I[i+1]:0) +
748  ((offd_i_I)?offd_i_I[i+1]:0);
749  offd_I[i + nrows + 1] = offd_I[i+1] + offd_r_nnz + offd_i_nnz;
750 
751  if (offd_r_I)
752  {
753  const int off_i = (offd_i_I)?(offd_i_I[i+1] - offd_i_I[i]):0;
754  for (int j=0; j<offd_r_I[i+1] - offd_r_I[i]; j++)
755  {
756  offd_J[offd_I[i] + j] =
757  cinvmap[cmapa[cmap_r[offd_r_J[offd_r_I[i] + j]]]];
758  offd_D[offd_I[i] + j] = offd_r_D[offd_r_I[i] + j];
759 
760  offd_J[offd_I[i+nrows] + off_i + j] =
761  cinvmap[cmapb[cmap_r[offd_r_J[offd_r_I[i] + j]]]];
762  offd_D[offd_I[i+nrows] + off_i + j] =
763  factor * offd_r_D[offd_r_I[i] + j];
764  }
765  }
766  if (offd_i_I)
767  {
768  const int off_r = (offd_r_I)?(offd_r_I[i+1] - offd_r_I[i]):0;
769  for (int j=0; j<offd_i_I[i+1] - offd_i_I[i]; j++)
770  {
771  offd_J[offd_I[i] + off_r + j] =
772  cinvmap[cmapb[cmap_i[offd_i_J[offd_i_I[i] + j]]]];
773  offd_D[offd_I[i] + off_r + j] = -offd_i_D[offd_i_I[i] + j];
774 
775  offd_J[offd_I[i+nrows] + j] =
776  cinvmap[cmapa[cmap_i[offd_i_J[offd_i_I[i] + j]]]];
777  offd_D[offd_I[i+nrows] + j] = factor * offd_i_D[offd_i_I[i] + j];
778  }
779  }
780  }
781 
782  // Construct the combined matrix
783  HypreParMatrix * A = new HypreParMatrix(comm_,
784  2 * global_num_rows,
785  2 * global_num_cols,
786  row_starts, col_starts,
787  diag_I, diag_J, diag_D,
788  offd_I, offd_J, offd_D,
789  2 * num_cols_offd, cmap,
790  true);
791 
792 #if MFEM_HYPRE_VERSION <= 22200
793  // Give the new matrix ownership of row_starts and col_starts
794  hypre_ParCSRMatrix *hA = (hypre_ParCSRMatrix*)(*A);
795 
796  hypre_ParCSRMatrixSetRowStartsOwner(hA,1);
797  hypre_ParCSRMatrixSetColStartsOwner(hA,1);
798 #else
799  mfem_hypre_TFree_host(row_starts);
800  mfem_hypre_TFree_host(col_starts);
801 #endif
802 
803  return A;
804 }
805 
806 void
807 ComplexHypreParMatrix::getColStartStop(const HypreParMatrix * A_r,
808  const HypreParMatrix * A_i,
809  int & num_recv_procs,
810  HYPRE_BigInt *& offd_col_start_stop
811  ) const
812 {
813  hypre_ParCSRCommPkg * comm_pkg_r =
814  (A_r) ? hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix*)(*A_r)) : NULL;
815  hypre_ParCSRCommPkg * comm_pkg_i =
816  (A_i) ? hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix*)(*A_i)) : NULL;
817 
818  std::set<HYPRE_Int> send_procs, recv_procs;
819  if ( comm_pkg_r )
820  {
821  for (HYPRE_Int i=0; i<comm_pkg_r->num_sends; i++)
822  {
823  send_procs.insert(comm_pkg_r->send_procs[i]);
824  }
825  for (HYPRE_Int i=0; i<comm_pkg_r->num_recvs; i++)
826  {
827  recv_procs.insert(comm_pkg_r->recv_procs[i]);
828  }
829  }
830  if ( comm_pkg_i )
831  {
832  for (HYPRE_Int i=0; i<comm_pkg_i->num_sends; i++)
833  {
834  send_procs.insert(comm_pkg_i->send_procs[i]);
835  }
836  for (HYPRE_Int i=0; i<comm_pkg_i->num_recvs; i++)
837  {
838  recv_procs.insert(comm_pkg_i->recv_procs[i]);
839  }
840  }
841 
842  num_recv_procs = (int)recv_procs.size();
843 
844  HYPRE_BigInt loc_start_stop[2];
845  offd_col_start_stop = new HYPRE_BigInt[2 * num_recv_procs];
846 
847  const HYPRE_BigInt * row_part = (A_r) ? A_r->RowPart() :
848  ((A_i) ? A_i->RowPart() : NULL);
849 
850  int row_part_ind = (HYPRE_AssumedPartitionCheck()) ? 0 : myid_;
851  loc_start_stop[0] = row_part[row_part_ind];
852  loc_start_stop[1] = row_part[row_part_ind+1];
853 
854  MPI_Request * req = new MPI_Request[send_procs.size()+recv_procs.size()];
855  MPI_Status * stat = new MPI_Status[send_procs.size()+recv_procs.size()];
856  int send_count = 0;
857  int recv_count = 0;
858  int tag = 0;
859 
860  std::set<HYPRE_Int>::iterator sit;
861  for (sit=send_procs.begin(); sit!=send_procs.end(); sit++)
862  {
863  MPI_Isend(loc_start_stop, 2, HYPRE_MPI_BIG_INT,
864  *sit, tag, comm_, &req[send_count]);
865  send_count++;
866  }
867  for (sit=recv_procs.begin(); sit!=recv_procs.end(); sit++)
868  {
869  MPI_Irecv(&offd_col_start_stop[2*recv_count], 2, HYPRE_MPI_BIG_INT,
870  *sit, tag, comm_, &req[send_count+recv_count]);
871  recv_count++;
872  }
873 
874  MPI_Waitall(send_count+recv_count, req, stat);
875 
876  delete [] req;
877  delete [] stat;
878 }
879 
880 #endif // MFEM_USE_MPI
881 
882 }
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:93
double Control[UMFPACK_CONTROL]
const double * HostReadData() const
Definition: sparsemat.hpp:277
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:1588
virtual Operator & real()
Real or imaginary part accessor methods.
virtual void Mult(const Vector &b, Vector &x) const
This is solving the system A x = b.
ComplexSparseMatrix * mat
ComplexHypreParMatrix(HypreParMatrix *A_Real, HypreParMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
Mimic the action of a complex operator using two real operators.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
virtual void MultTranspose(const Vector &b, Vector &x) const
This is solving the system: A^H x = b (when transa = false) This is equivalent to solving the transpo...
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Alternate convention for damping operators.
ComplexOperator(Operator *Op_Real, Operator *Op_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
Constructs complex operator object.
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:457
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
HypreParMatrix * GetSystemMatrix() const
Data type sparse matrix.
Definition: sparsemat.hpp:50
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double b
Definition: lissajous.cpp:42
const int * HostReadJ() const
Definition: sparsemat.hpp:261
SparseMatrix * GetSystemMatrix() const
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1481
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
virtual SparseMatrix & real()
Real or imaginary part accessor methods.
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:534
Specialization of the ComplexOperator built from a pair of Sparse Matrices.
HYPRE_Int HYPRE_BigInt
virtual Operator & imag()
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:154
const int * HostReadI() const
Definition: sparsemat.hpp:245
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a ComplexSparseMatrix.
Convention GetConvention() const
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
Definition: vector.hpp:238
HYPRE_BigInt * ColPart()
Returns the column partitioning.
Definition: hypre.hpp:573
HYPRE_BigInt * RowPart()
Returns the row partitioning.
Definition: hypre.hpp:569
virtual SparseMatrix & imag()
Native convention for Hermitian operators.
Vector data type.
Definition: vector.hpp:58
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:581
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
HYPRE_BigInt GetGlobalNumCols() const
Return the global number of columns.
Definition: hypre.hpp:639
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void GetOffd(SparseMatrix &offd, HYPRE_BigInt *&cmap) const
Get the local off-diagonal block. NOTE: &#39;offd&#39; will not own any data.
Definition: hypre.cpp:1557
virtual HypreParMatrix & imag()
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
virtual HypreParMatrix & real()
Real or imaginary part accessor methods.
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:473
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition: hypre.hpp:635
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28