MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
superlu.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 
14 #ifdef MFEM_USE_SUPERLU
15 #ifdef MFEM_USE_MPI
16 
17 #include "superlu.hpp"
18 
19 // SuperLU headers
20 #include "superlu_defs.h"
21 #include "superlu_ddefs.h"
22 
23 #if XSDK_INDEX_SIZE == 64
24 #error "SuperLUDist has been built with 64bit integers. This is not supported"
25 #endif
26 
27 // For now, it is assumed that HYPRE_BigInt is int.
28 #if defined(HYPRE_BIGINT) || defined(HYPRE_MIXEDINT)
29 #error "SuperLUDist support requires HYPRE_BigInt == int, for now."
30 #endif
31 
32 #if SUPERLU_DIST_MAJOR_VERSION > 6 || \
33  (SUPERLU_DIST_MAJOR_VERSION == 6 && SUPERLU_DIST_MINOR_VERSION > 2)
34 #define ScalePermstruct_t dScalePermstruct_t
35 #define LUstruct_t dLUstruct_t
36 #define SOLVEstruct_t dSOLVEstruct_t
37 #define ScalePermstructFree dScalePermstructFree
38 #define Destroy_LU dDestroy_LU
39 #define LUstructFree dLUstructFree
40 #define LUstructInit dLUstructInit
41 #endif
42 
43 
44 using namespace std;
45 
46 namespace mfem
47 {
48 unsigned int superlu_internal::sqrti( const unsigned int & a )
49 {
50  unsigned int a_ = a;
51  unsigned int rem = 0;
52  unsigned int root = 0;
53  unsigned short len = sizeof(int); len <<= 2;
54  unsigned short shift = (unsigned short)((len<<1) - 2);
55 
56  for (int i=0; i<len; i++)
57  {
58  root <<= 1;
59  rem = ((rem << 2) + (a_ >> shift));
60  a_ <<= 2;
61  root ++;
62  if (root <= rem)
63  {
64  rem -= root;
65  root++;
66  }
67  else
68  {
69  root--;
70  }
71  }
72  return (root >> 1);
73 }
74 
75 SuperLURowLocMatrix::SuperLURowLocMatrix(MPI_Comm comm,
76  int num_loc_rows, int first_loc_row,
77  int glob_nrows, int glob_ncols,
78  int *I, int *J, double *data)
79  : comm_(comm),
80  rowLocPtr_(NULL)
81 {
82  // Set mfem::Operator member data
83  height = num_loc_rows;
84  width = num_loc_rows;
85 
86  // Allocate SuperLU's SuperMatrix struct
87  rowLocPtr_ = new SuperMatrix;
88  SuperMatrix * A = (SuperMatrix*)rowLocPtr_;
89 
90  A->Store = NULL;
91 
92  int m = glob_nrows;
93  int n = glob_ncols;
94  int nnz_loc = I[num_loc_rows];
95  int m_loc = num_loc_rows;
96  int fst_row = first_loc_row;
97 
98  double * nzval = NULL;
99  int * colind = NULL;
100  int * rowptr = NULL;
101 
102  if ( !(nzval = doubleMalloc_dist(nnz_loc)) )
103  {
104  ABORT("Malloc fails for nzval[].");
105  }
106  for (int i=0; i<nnz_loc; i++)
107  {
108  nzval[i] = data[i];
109  }
110 
111  if ( !(colind = intMalloc_dist(nnz_loc)) )
112  {
113  ABORT("Malloc fails for colind[].");
114  }
115  for (int i=0; i<nnz_loc; i++)
116  {
117  colind[i] = J[i];
118  }
119 
120  if ( !(rowptr = intMalloc_dist(m_loc+1)) )
121  {
122  ABORT("Malloc fails for rowptr[].");
123  }
124  for (int i=0; i<=m_loc; i++)
125  {
126  rowptr[i] = I[i];
127  }
128 
129  // Assign he matrix data to SuperLU's SuperMatrix structure
130  dCreate_CompRowLoc_Matrix_dist(A, m, n, nnz_loc, m_loc, fst_row,
131  nzval, colind, rowptr,
132  SLU_NR_loc, SLU_D, SLU_GE);
133 }
134 
136  : comm_(hypParMat.GetComm()),
137  rowLocPtr_(NULL)
138 {
139  rowLocPtr_ = new SuperMatrix;
140  SuperMatrix * A = (SuperMatrix*)rowLocPtr_;
141 
142  A->Store = NULL;
143 
144  // First cast the parameter to a hypre_ParCSRMatrix
145  hypre_ParCSRMatrix * parcsr_op =
146  (hypre_ParCSRMatrix *)const_cast<HypreParMatrix&>(hypParMat);
147 
148  MFEM_ASSERT(parcsr_op != NULL,"SuperLU: const_cast failed in SetOperator");
149 
150  // Create the SuperMatrix A by borrowing the internal data from a
151  // hypre_CSRMatrix.
152  hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
153  hypre_CSRMatrixSetDataOwner(csr_op,0);
154 #if MFEM_HYPRE_VERSION >= 21600
155  // For now, this method assumes that HYPRE_BigInt is int. Also, csr_op->num_cols
156  // is of type HYPRE_Int, so if we want to check for big indices in
157  // csr_op->big_j, we'll have to check all entries and that check will only be
158  // necessary in HYPRE_MIXEDINT mode which is not supported at the moment.
159  hypre_CSRMatrixBigJtoJ(csr_op);
160 #endif
161 
162  int m = parcsr_op->global_num_rows;
163  int n = parcsr_op->global_num_cols;
164  int fst_row = parcsr_op->first_row_index;
165  int nnz_loc = csr_op->num_nonzeros;
166  int m_loc = csr_op->num_rows;
167 
168  height = m_loc;
169  width = m_loc;
170 
171  double * nzval = csr_op->data;
172  int * colind = csr_op->j;
173  int * rowptr = NULL;
174 
175  // The "i" array cannot be stolen from the hypre_CSRMatrix so we'll copy it
176  if ( !(rowptr = intMalloc_dist(m_loc+1)) )
177  {
178  ABORT("Malloc fails for rowptr[].");
179  }
180  for (int i=0; i<=m_loc; i++)
181  {
182  rowptr[i] = (csr_op->i)[i];
183  }
184 
185  // Everything has been copied or abducted so delete the structure
186  hypre_CSRMatrixDestroy(csr_op);
187 
188  // Assign he matrix data to SuperLU's SuperMatrix structure
189  dCreate_CompRowLoc_Matrix_dist(A, m, n, nnz_loc, m_loc, fst_row,
190  nzval, colind, rowptr,
191  SLU_NR_loc, SLU_D, SLU_GE);
192 
193  // Save global number of columns (width) of the matrix
194  num_global_cols = n;
195 }
196 
198 {
199  SuperMatrix * A = (SuperMatrix*)rowLocPtr_;
200 
201  // Delete the internal data
202  Destroy_CompRowLoc_Matrix_dist(A);
203 
204  // Delete the struct
205  if ( A != NULL ) { delete A; }
206 }
207 
209  : comm_(comm),
210  APtr_(NULL),
211  optionsPtr_(NULL),
212  statPtr_(NULL),
213  ScalePermstructPtr_(NULL),
214  LUstructPtr_(NULL),
215  SOLVEstructPtr_(NULL),
216  gridPtr_(NULL),
217  berr_(NULL),
218  perm_r_(NULL),
219  nrhs_(1),
220  nprow_(0),
221  npcol_(0),
222  firstSolveWithThisA_(false),
223  gridInitialized_(false),
224  LUStructInitialized_(false)
225 {
226  this->Init();
227 }
228 
230  : comm_(A.GetComm()),
231  APtr_(&A),
232  optionsPtr_(NULL),
233  statPtr_(NULL),
234  ScalePermstructPtr_(NULL),
235  LUstructPtr_(NULL),
236  SOLVEstructPtr_(NULL),
237  gridPtr_(NULL),
238  berr_(NULL),
239  perm_r_(NULL),
240  nrhs_(1),
241  nprow_(0),
242  npcol_(0),
243  firstSolveWithThisA_(true),
244  gridInitialized_(false),
245  LUStructInitialized_(false)
246 {
247  height = A.Height();
248  width = A.Width();
249 
250  this->Init();
251 }
252 
254 {
255  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
256  SuperLUStat_t * stat = (SuperLUStat_t*)statPtr_;
257  ScalePermstruct_t * SPstruct = (ScalePermstruct_t*)ScalePermstructPtr_;
258  LUstruct_t * LUstruct = (LUstruct_t*)LUstructPtr_;
259  SOLVEstruct_t * SOLVEstruct = (SOLVEstruct_t*)SOLVEstructPtr_;
260  gridinfo_t * grid = (gridinfo_t*)gridPtr_;
261 
262  SUPERLU_FREE(berr_);
263  PStatFree(stat);
264 
265  if ( LUStructInitialized_ )
266  {
267  ScalePermstructFree(SPstruct);
268  Destroy_LU(APtr_->GetGlobalNumColumns(), grid, LUstruct);
269  LUstructFree(LUstruct);
270  }
271 
272  if ( options->SolveInitialized )
273  {
274  dSolveFinalize(options, SOLVEstruct);
275  }
276 
277  if ( options != NULL ) { delete options; }
278  if ( stat != NULL ) { delete stat; }
279  if ( SPstruct != NULL ) { delete SPstruct; }
280  if ( LUstruct != NULL ) { delete LUstruct; }
281  if ( SOLVEstruct != NULL ) { delete SOLVEstruct; }
282  if ( grid != NULL ) { delete grid; }
283  if ( perm_r_ != NULL ) { SUPERLU_FREE(perm_r_); }
284 }
285 
286 void SuperLUSolver::Init()
287 {
288  MPI_Comm_size(comm_, &numProcs_);
289  MPI_Comm_rank(comm_, &myid_);
290 
291  optionsPtr_ = new superlu_dist_options_t;
292  statPtr_ = new SuperLUStat_t;
293  ScalePermstructPtr_ = new ScalePermstruct_t;
294  LUstructPtr_ = new LUstruct_t;
295  SOLVEstructPtr_ = new SOLVEstruct_t;
296  gridPtr_ = new gridinfo_t;
297 
298  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
299  SuperLUStat_t * stat = (SuperLUStat_t*)statPtr_;
300 
301  if ( !(berr_ = doubleMalloc_dist(nrhs_)) )
302  {
303  ABORT("Malloc fails for berr[].");
304  }
305 
306  // Set default options
307  set_default_options_dist(options);
308 
309  // Choose nprow and npcol so that the process grid is as square as possible.
310  // If the processes cannot be divided evenly, keep the row dimension smaller
311  // than the column dimension.
312 
313  nprow_ = (int)superlu_internal::sqrti((unsigned int)numProcs_);
314  while (numProcs_ % nprow_ != 0 && nprow_ > 0)
315  {
316  nprow_--;
317  }
318 
319  npcol_ = (int)(numProcs_ / nprow_);
320  MFEM_ASSERT(nprow_ * npcol_ == numProcs_, "");
321 
322  PStatInit(stat); // Initialize the statistics variables.
323 }
324 
325 void SuperLUSolver::SetPrintStatistics( bool print_stat )
326 {
327  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
328 
329  yes_no_t opt = print_stat?YES:NO;
330 
331  options->PrintStat = opt;
332 }
333 
335 {
336  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
337 
338  yes_no_t opt = equil?YES:NO;
339 
340  options->Equil = opt;
341 }
342 
344 {
345  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
346 
347  colperm_t opt = (colperm_t)col_perm;
348 
349  options->ColPerm = opt;
350 }
351 
353  Array<int> * perm )
354 {
355  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
356 
357  rowperm_t opt = (rowperm_t)row_perm;
358 
359  options->RowPerm = opt;
360 
361  if ( opt == MY_PERMR )
362  {
363  if ( perm == NULL )
364  {
365  mfem_error("SuperLUSolver::SetRowPermutation :"
366  " permutation vector not set!");
367  }
368 
369  if ( !(perm_r_ = intMalloc_dist(perm->Size())) )
370  {
371  ABORT("Malloc fails for perm_r[].");
372  }
373  for (int i=0; i<perm->Size(); i++)
374  {
375  perm_r_[i] = (*perm)[i];
376  }
377  }
378 }
379 
381 {
382  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
383 
384  trans_t opt = (trans_t)trans;
385 
386  options->Trans = opt;
387 }
388 
390 {
391  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
392 
393  IterRefine_t opt = (IterRefine_t)iter_ref;
394 
395  options->IterRefine = opt;
396 }
397 
399 {
400  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
401 
402  yes_no_t opt = rtp?YES:NO;
403 
404  options->ReplaceTinyPivot = opt;
405 }
406 
407 void SuperLUSolver::SetNumLookAheads( int num_lookaheads )
408 {
409  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
410 
411  options->num_lookaheads = num_lookaheads;
412 }
413 
415 {
416  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
417 
418  yes_no_t opt = etree?YES:NO;
419 
420  options->lookahead_etree = opt;
421 }
422 
424 {
425  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
426 
427  yes_no_t opt = sym?YES:NO;
428 
429  options->SymPattern = opt;
430 }
431 
433 {
434  gridinfo_t * grid = (gridinfo_t*)gridPtr_;
435 
436  // Make sure the values of nprow and npcol are reasonable
437  if ( ((nprow_ * npcol_) > numProcs_) || ((nprow_ * npcol_) < 1) )
438  {
439  if ( myid_ == 0 )
440  {
441  mfem::err << "Warning: User specified nprow and npcol are such that "
442  << "(nprow * npcol) > numProcs or (nprow * npcol) < 1. "
443  << "Using default values for nprow and npcol instead."
444  << endl;
445  }
446 
447  nprow_ = (int)superlu_internal::sqrti((unsigned int)numProcs_);
448  while (numProcs_ % nprow_ != 0 && nprow_ > 0)
449  {
450  nprow_--;
451  }
452 
453  npcol_ = (int)(numProcs_ / nprow_);
454  MFEM_ASSERT(nprow_ * npcol_ == numProcs_, "");
455  }
456 
457  superlu_gridinit(comm_, nprow_, npcol_, grid);
458 
459  gridInitialized_ = true;
460 }
461 
463 {
464  if ( gridInitialized_ )
465  {
466  gridinfo_t * grid = (gridinfo_t*)gridPtr_;
467 
468  superlu_gridexit(grid);
469  }
470 
471  gridInitialized_ = false;
472 }
473 
474 void SuperLUSolver::Mult( const Vector & x, Vector & y ) const
475 {
476  MFEM_ASSERT(APtr_ != NULL,
477  "SuperLU Error: The operator must be set before"
478  " the system can be solved.");
479 
480  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
481  SuperLUStat_t * stat = (SuperLUStat_t*)statPtr_;
482  SuperMatrix * A = (SuperMatrix*)APtr_->InternalData();
483 
484  ScalePermstruct_t * SPstruct = (ScalePermstruct_t*)ScalePermstructPtr_;
485  LUstruct_t * LUstruct = (LUstruct_t*)LUstructPtr_;
486  SOLVEstruct_t * SOLVEstruct = (SOLVEstruct_t*)SOLVEstructPtr_;
487  gridinfo_t * grid = (gridinfo_t*)gridPtr_;
488 
490  {
491  options->Fact = FACTORED; // Indicate the factored form of A is supplied.
492  }
493  else // This is the first solve with this A
494  {
495  firstSolveWithThisA_ = false;
496 
497  // Make sure that the parameters have been initialized The only parameter
498  // we might have to worry about is ScalePermstruct, if the user is
499  // supplying a row or column permutation.
500 
501  // Initialize ScalePermstruct and LUstruct.
502  SPstruct->DiagScale = NOEQUIL;
503 
504  // Transfer ownership of the row permutations if available
505  if ( perm_r_ != NULL )
506  {
507  SPstruct->perm_r = perm_r_;
508  perm_r_ = NULL;
509  }
510  else
511  {
512  if ( !(SPstruct->perm_r = intMalloc_dist(A->nrow)) )
513  {
514  ABORT("Malloc fails for perm_r[].");
515  }
516  }
517  if ( !(SPstruct->perm_c = intMalloc_dist(A->ncol)) )
518  {
519  ABORT("Malloc fails for perm_c[].");
520  }
521 
522  LUstructInit(A->ncol, LUstruct);
523  LUStructInitialized_ = true;
524  }
525 
526  // SuperLU overwrites x with y, so copy x to y and pass that to the solve
527  // routine.
528 
529  const double *xPtr = x.HostRead();
530  y = xPtr;
531  double * yPtr = y.HostReadWrite();
532 
533  int info = -1, locSize = y.Size();
534 
535  // Solve the system
536  pdgssvx(options, A, SPstruct, yPtr, locSize, nrhs_, grid,
537  LUstruct, SOLVEstruct, berr_, stat, &info);
538 
539  if ( info != 0 )
540  {
541  if ( info < 0 )
542  {
543  switch (-info)
544  {
545  case 1:
546  MFEM_ABORT("SuperLU: SuperLU options are invalid.");
547  break;
548  case 2:
549  MFEM_ABORT("SuperLU: Matrix A (in Ax=b) is invalid.");
550  break;
551  case 5:
552  MFEM_ABORT("SuperLU: Vector b dimension (in Ax=b) is invalid.");
553  break;
554  case 6:
555  MFEM_ABORT("SuperLU: Number of right-hand sides is invalid.");
556  break;
557  default:
558  MFEM_ABORT("SuperLU: Parameter with index "
559  << -info << "invalid. (1-indexed)");
560  break;
561  }
562  }
563  else if ( info <= A->ncol )
564  {
565  MFEM_ABORT("SuperLU: Found a singular matrix, U("
566  << info << "," << info << ") is exactly zero.");
567  }
568  else if ( info > A->ncol )
569  {
570  MFEM_ABORT("SuperLU: Memory allocation error with "
571  << info - A->ncol << " bytes already allocated,");
572  }
573  else
574  {
575  MFEM_ABORT("Unknown SuperLU Error");
576  }
577  }
578 }
579 
581 {
582  // Verify that we have a compatible operator
583  APtr_ = dynamic_cast<const SuperLURowLocMatrix*>(&op);
584  if ( APtr_ == NULL )
585  {
586  mfem_error("SuperLUSolver::SetOperator : not SuperLURowLocMatrix!");
587  }
588 
589  // Everything is OK so finish setting the operator
590  firstSolveWithThisA_ = true;
591 
592  // Set mfem::Operator member data
593  height = op.Height();
594  width = op.Width();
595 
596  // Initialize the processor grid if necessary
597  if (!gridInitialized_)
598  {
599  this->SetupGrid();
600  }
601 }
602 
603 } // mfem namespace
604 
605 #endif // MFEM_USE_MPI
606 #endif // MFEM_USE_SUPERLU
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void * InternalData() const
Definition: superlu.hpp:73
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: superlu.cpp:474
void * ScalePermstructPtr_
Definition: superlu.hpp:147
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void SetRowPermutation(superlu::RowPerm row_perm, Array< int > *perm=NULL)
Definition: superlu.cpp:352
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
void SetSymmetricPattern(bool sym)
Definition: superlu.cpp:423
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: superlu.cpp:580
HYPRE_BigInt GetGlobalNumColumns() const
Definition: superlu.hpp:75
void SetIterativeRefine(superlu::IterRefine iter_ref)
Definition: superlu.cpp:389
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void SetColumnPermutation(superlu::ColPerm col_perm)
Definition: superlu.cpp:343
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
void SetNumLookAheads(int num_lookaheads)
Definition: superlu.cpp:407
void SetPrintStatistics(bool print_stat)
Definition: superlu.cpp:325
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:442
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
A class to initialize the size of a Tensor.
Definition: dtensor.hpp:54
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
double a
Definition: lissajous.cpp:41
void * SOLVEstructPtr_
Definition: superlu.hpp:149
void SetLookAheadElimTree(bool etree)
Definition: superlu.cpp:414
unsigned int sqrti(const unsigned int &a)
Definition: superlu.cpp:48
void SetReplaceTinyPivot(bool rtp)
Definition: superlu.cpp:398
Vector data type.
Definition: vector.hpp:60
void SetTranspose(superlu::Trans trans)
Definition: superlu.cpp:380
const SuperLURowLocMatrix * APtr_
Definition: superlu.hpp:137
SuperLURowLocMatrix(MPI_Comm comm, int num_loc_rows, int first_loc_row, int glob_nrows, int glob_ncols, int *I, int *J, double *data)
Definition: superlu.cpp:75
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
SuperLUSolver(MPI_Comm comm)
Definition: superlu.cpp:208
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:458
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void SetEquilibriate(bool equil)
Definition: superlu.cpp:334