MFEM  v3.3.2
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../config/config.hpp"
13 
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 using namespace std;
24 
25 namespace mfem
26 {
27 unsigned int superlu_internal::sqrti( const unsigned int & a )
28 {
29  unsigned int a_ = a;
30  unsigned int rem = 0;
31  unsigned int root = 0;
32  unsigned short len = sizeof(int); len <<= 2;
33  unsigned short shift = (unsigned short)((len<<1) - 2);
34 
35  for (int i=0; i<len; i++)
36  {
37  root <<= 1;
38  rem = ((rem << 2) + (a_ >> shift));
39  a_ <<= 2;
40  root ++;
41  if (root <= rem)
42  {
43  rem -= root;
44  root++;
45  }
46  else
47  {
48  root--;
49  }
50  }
51  return (root >> 1);
52 }
53 
54 SuperLURowLocMatrix::SuperLURowLocMatrix(MPI_Comm comm,
55  int num_loc_rows, int first_loc_row,
56  int glob_nrows, int glob_ncols,
57  int *I, int *J, double *data)
58  : comm_(comm),
59  rowLocPtr_(NULL)
60 {
61  // Set mfem::Operator member data
62  height = num_loc_rows;
63  width = num_loc_rows;
64 
65  // Allocate SuperLU's SuperMatrix struct
66  rowLocPtr_ = new SuperMatrix;
67  SuperMatrix * A = (SuperMatrix*)rowLocPtr_;
68 
69  A->Store = NULL;
70 
71  int m = glob_nrows;
72  int n = glob_ncols;
73  int nnz_loc = I[num_loc_rows];
74  int m_loc = num_loc_rows;
75  int fst_row = first_loc_row;
76 
77  double * nzval = NULL;
78  int * colind = NULL;
79  int * rowptr = NULL;
80 
81  if ( !(nzval = doubleMalloc_dist(nnz_loc)) )
82  {
83  ABORT("Malloc fails for nzval[].");
84  }
85  for (int i=0; i<nnz_loc; i++)
86  {
87  nzval[i] = data[i];
88  }
89 
90  if ( !(colind = intMalloc_dist(nnz_loc)) )
91  {
92  ABORT("Malloc fails for colind[].");
93  }
94  for (int i=0; i<nnz_loc; i++)
95  {
96  colind[i] = J[i];
97  }
98 
99  if ( !(rowptr = intMalloc_dist(m_loc+1)) )
100  {
101  ABORT("Malloc fails for rowptr[].");
102  }
103  for (int i=0; i<=m_loc; i++)
104  {
105  rowptr[i] = I[i];
106  }
107 
108  // Assign he matrix data to SuperLU's SuperMatrix structure
109  dCreate_CompRowLoc_Matrix_dist(A, m, n, nnz_loc, m_loc, fst_row,
110  nzval, colind, rowptr,
111  SLU_NR_loc, SLU_D, SLU_GE);
112 }
113 
115  : comm_(hypParMat.GetComm()),
116  rowLocPtr_(NULL)
117 {
118  rowLocPtr_ = new SuperMatrix;
119  SuperMatrix * A = (SuperMatrix*)rowLocPtr_;
120 
121  A->Store = NULL;
122 
123  // First cast the parameter to a hypre_ParCSRMatrix
124  hypre_ParCSRMatrix * parcsr_op =
125  (hypre_ParCSRMatrix *)const_cast<HypreParMatrix&>(hypParMat);
126 
127  MFEM_ASSERT(parcsr_op != NULL,"SuperLU: const_cast failed in SetOperator");
128 
129  // Create the SuperMatrix A by borrowing the internal data from a
130  // hypre_CSRMatrix.
131  hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
132  hypre_CSRMatrixSetDataOwner(csr_op,0);
133 
134  int m = parcsr_op->global_num_rows;
135  int n = parcsr_op->global_num_cols;
136  int fst_row = parcsr_op->first_row_index;
137  int nnz_loc = csr_op->num_nonzeros;
138  int m_loc = csr_op->num_rows;
139 
140  height = m_loc;
141  width = m_loc;
142 
143  double * nzval = csr_op->data;
144  int * colind = csr_op->j;
145  int * rowptr = NULL;
146 
147  // The "i" array cannot be stolen from the hypre_CSRMatrix so we'll copy it
148  if ( !(rowptr = intMalloc_dist(m_loc+1)) )
149  {
150  ABORT("Malloc fails for rowptr[].");
151  }
152  for (int i=0; i<=m_loc; i++)
153  {
154  rowptr[i] = (csr_op->i)[i];
155  }
156 
157  // Everything has been copied or abducted so delete the structure
158  hypre_CSRMatrixDestroy(csr_op);
159 
160  // Assign he matrix data to SuperLU's SuperMatrix structure
161  dCreate_CompRowLoc_Matrix_dist(A, m, n, nnz_loc, m_loc, fst_row,
162  nzval, colind, rowptr,
163  SLU_NR_loc, SLU_D, SLU_GE);
164 }
165 
167 {
168  SuperMatrix * A = (SuperMatrix*)rowLocPtr_;
169 
170  // Delete the internal data
171  Destroy_CompRowLoc_Matrix_dist(A);
172 
173  // Delete the struct
174  if ( A != NULL ) { delete A; }
175 }
176 
178  : comm_(comm),
179  APtr_(NULL),
180  optionsPtr_(NULL),
181  statPtr_(NULL),
182  ScalePermstructPtr_(NULL),
183  LUstructPtr_(NULL),
184  SOLVEstructPtr_(NULL),
185  gridPtr_(NULL),
186  berr_(NULL),
187  perm_r_(NULL),
188  nrhs_(1),
189  nprow_(0),
190  npcol_(0),
191  firstSolveWithThisA_(false),
192  gridInitialized_(false),
193  LUStructInitialized_(false)
194 {
195  this->Init();
196 }
197 
199  : comm_(A.GetComm()),
200  APtr_(&A),
201  optionsPtr_(NULL),
202  statPtr_(NULL),
203  ScalePermstructPtr_(NULL),
204  LUstructPtr_(NULL),
205  SOLVEstructPtr_(NULL),
206  gridPtr_(NULL),
207  berr_(NULL),
208  perm_r_(NULL),
209  nrhs_(1),
210  nprow_(0),
211  npcol_(0),
212  firstSolveWithThisA_(true),
213  gridInitialized_(false),
214  LUStructInitialized_(false)
215 {
216  height = A.Height();
217  width = A.Width();
218 
219  this->Init();
220 }
221 
223 {
224  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
225  SuperLUStat_t * stat = (SuperLUStat_t*)statPtr_;
226  ScalePermstruct_t * SPstruct = (ScalePermstruct_t*)ScalePermstructPtr_;
227  LUstruct_t * LUstruct = (LUstruct_t*)LUstructPtr_;
228  SOLVEstruct_t * SOLVEstruct = (SOLVEstruct_t*)SOLVEstructPtr_;
229  gridinfo_t * grid = (gridinfo_t*)gridPtr_;
230 
231  SUPERLU_FREE(berr_);
232  PStatFree(stat);
233 
234  if ( LUStructInitialized_ )
235  {
236  ScalePermstructFree(SPstruct);
237  Destroy_LU(width, grid, LUstruct);
238  LUstructFree(LUstruct);
239  }
240 
241  if ( options->SolveInitialized )
242  {
243  dSolveFinalize(options, SOLVEstruct);
244  }
245 
246  if ( options != NULL ) { delete options; }
247  if ( stat != NULL ) { delete stat; }
248  if ( SPstruct != NULL ) { delete SPstruct; }
249  if ( LUstruct != NULL ) { delete LUstruct; }
250  if ( SOLVEstruct != NULL ) { delete SOLVEstruct; }
251  if ( grid != NULL ) { delete grid; }
252  if ( perm_r_ != NULL ) { SUPERLU_FREE(perm_r_); }
253 }
254 
255 void SuperLUSolver::Init()
256 {
257  MPI_Comm_size(comm_, &numProcs_);
258  MPI_Comm_rank(comm_, &myid_);
259 
260  optionsPtr_ = new superlu_dist_options_t;
261  statPtr_ = new SuperLUStat_t;
262  ScalePermstructPtr_ = new ScalePermstruct_t;
263  LUstructPtr_ = new LUstruct_t;
264  SOLVEstructPtr_ = new SOLVEstruct_t;
265  gridPtr_ = new gridinfo_t;
266 
267  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
268  SuperLUStat_t * stat = (SuperLUStat_t*)statPtr_;
269 
270  if ( !(berr_ = doubleMalloc_dist(nrhs_)) )
271  {
272  ABORT("Malloc fails for berr[].");
273  }
274 
275  // Set default options
276  set_default_options_dist(options);
277 
278  options->ParSymbFact = YES;
279  options->ColPerm = NATURAL;
280 
281  // Choose nprow and npcol so that the process grid is as square as possible.
282  // If the processes cannot be divided evenly, keep the row dimension smaller
283  // than the column dimension.
284 
285  nprow_ = (int)superlu_internal::sqrti((unsigned int)numProcs_);
286  while (numProcs_ % nprow_ != 0 && nprow_ > 0)
287  {
288  nprow_--;
289  }
290 
291  npcol_ = (int)(numProcs_ / nprow_);
292  assert(nprow_ * npcol_ == numProcs_);
293 
294  PStatInit(stat); // Initialize the statistics variables.
295 }
296 
297 void SuperLUSolver::SetPrintStatistics( bool print_stat )
298 {
299  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
300 
301  yes_no_t opt = print_stat?YES:NO;
302 
303  options->PrintStat = opt;
304 }
305 
307 {
308  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
309 
310  yes_no_t opt = equil?YES:NO;
311 
312  options->Equil = opt;
313 }
314 
316 {
317  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
318 
319  colperm_t opt = (colperm_t)col_perm;
320 
321  options->ColPerm = opt;
322 }
323 
325  Array<int> * perm )
326 {
327  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
328 
329  rowperm_t opt = (rowperm_t)row_perm;
330 
331  options->RowPerm = opt;
332 
333  if ( opt == MY_PERMR )
334  {
335  if ( perm == NULL )
336  {
337  mfem_error("SuperLUSolver::SetRowPermutation :"
338  " permutation vector not set!");
339  }
340 
341  if ( !(perm_r_ = intMalloc_dist(perm->Size())) )
342  {
343  ABORT("Malloc fails for perm_r[].");
344  }
345  for (int i=0; i<perm->Size(); i++)
346  {
347  perm_r_[i] = (*perm)[i];
348  }
349  }
350 }
351 
353 {
354  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
355 
356  trans_t opt = (trans_t)trans;
357 
358  options->Trans = opt;
359 }
360 
362 {
363  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
364 
365  IterRefine_t opt = (IterRefine_t)iter_ref;
366 
367  options->IterRefine = opt;
368 }
369 
371 {
372  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
373 
374  yes_no_t opt = rtp?YES:NO;
375 
376  options->ReplaceTinyPivot = opt;
377 }
378 
379 void SuperLUSolver::SetNumLookAheads( int num_lookaheads )
380 {
381  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
382 
383  options->num_lookaheads = num_lookaheads;
384 }
385 
387 {
388  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
389 
390  yes_no_t opt = etree?YES:NO;
391 
392  options->lookahead_etree = opt;
393 }
394 
396 {
397  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
398 
399  yes_no_t opt = sym?YES:NO;
400 
401  options->SymPattern = opt;
402 }
403 
405 {
406  gridinfo_t * grid = (gridinfo_t*)gridPtr_;
407 
408  // Make sure the values of nprow and npcol are reasonable
409  if ( ((nprow_ * npcol_) > numProcs_) || ((nprow_ * npcol_) < 1) )
410  {
411  if ( myid_ == 0 )
412  {
413  mfem:err << "Warning: User specified nprow and npcol are such that "
414  << "(nprow * npcol) > numProcs or (nprow * npcol) < 1. "
415  << "Using default values for nprow and npcol instead." << endl;
416  }
417 
418  nprow_ = (int)superlu_internal::sqrti((unsigned int)numProcs_);
419  while (numProcs_ % nprow_ != 0 && nprow_ > 0)
420  {
421  nprow_--;
422  }
423 
424  npcol_ = (int)(numProcs_ / nprow_);
425  assert(nprow_ * npcol_ == numProcs_);
426  }
427 
428  superlu_gridinit(comm_, nprow_, npcol_, grid);
429 
430  gridInitialized_ = true;
431 }
432 
434 {
435  if ( gridInitialized_ )
436  {
437  gridinfo_t * grid = (gridinfo_t*)gridPtr_;
438 
439  superlu_gridexit(grid);
440  }
441 
442  gridInitialized_ = false;
443 }
444 
445 void SuperLUSolver::Mult( const Vector & x, Vector & y ) const
446 {
447  MFEM_ASSERT(APtr_ != NULL,
448  "SuperLU Error: The operator must be set before"
449  " the system can be solved.");
450 
451  superlu_dist_options_t * options = (superlu_dist_options_t*)optionsPtr_;
452  SuperLUStat_t * stat = (SuperLUStat_t*)statPtr_;
453  SuperMatrix * A = (SuperMatrix*)APtr_->InternalData();
454 
455  ScalePermstruct_t * SPstruct = (ScalePermstruct_t*)ScalePermstructPtr_;
456  LUstruct_t * LUstruct = (LUstruct_t*)LUstructPtr_;
457  SOLVEstruct_t * SOLVEstruct = (SOLVEstruct_t*)SOLVEstructPtr_;
458  gridinfo_t * grid = (gridinfo_t*)gridPtr_;
459 
461  {
462  options->Fact = FACTORED; // Indicate the factored form of A is supplied.
463  }
464  else // This is the first solve with this A
465  {
466  firstSolveWithThisA_ = false;
467 
468  // Make sure that the parameters have been initialized The only parameter
469  // we might have to worry about is ScalePermstruct, if the user is
470  // supplying a row or column permutation.
471 
472  // Initialize ScalePermstruct and LUstruct.
473  SPstruct->DiagScale = NOEQUIL;
474 
475  // Transfer ownership of the row permutations if available
476  if ( perm_r_ != NULL )
477  {
478  SPstruct->perm_r = perm_r_;
479  perm_r_ = NULL;
480  }
481  else
482  {
483  if ( !(SPstruct->perm_r = intMalloc_dist(A->nrow)) )
484  {
485  ABORT("Malloc fails for perm_r[].");
486  }
487  }
488  if ( !(SPstruct->perm_c = intMalloc_dist(A->ncol)) )
489  {
490  ABORT("Malloc fails for perm_c[].");
491  }
492 
493  LUstructInit(A->ncol, LUstruct);
494  LUStructInitialized_ = true;
495  }
496 
497  // SuperLU overwrites x with y, so copy x to y and pass that to the solve
498  // routine.
499 
500  y = x;
501 
502  double* yPtr = (double*)y;
503  int info = -1, locSize = y.Size();
504 
505  // Solve the system
506  pdgssvx(options, A, SPstruct, yPtr, locSize, nrhs_, grid,
507  LUstruct, SOLVEstruct, berr_, stat, &info);
508 
509  if ( info != 0 )
510  {
511  if ( info <= A->ncol )
512  {
513  MFEM_ABORT("SuperLU: Found a singular matrix, U("
514  << info << "," << info << ") is exactly zero.");
515  }
516  else if ( info > A->ncol )
517  {
518  MFEM_ABORT("SuperLU: Memory allocation error with "
519  << info - A->ncol << " bytes already allocated,");
520  }
521  else
522  {
523  MFEM_ABORT("Unknown SuperLU Error");
524  }
525  }
526 }
527 
529 {
530  // Verify that we have a compatible operator
531  APtr_ = dynamic_cast<const SuperLURowLocMatrix*>(&op);
532  if ( APtr_ == NULL )
533  {
534  mfem_error("SuperLUSolver::SetOperator : not SuperLURowLocMatrix!");
535  }
536 
537  // Everything is OK so finish setting the operator
538  firstSolveWithThisA_ = true;
539 
540  // Set mfem::Operator member data
541  height = op.Height();
542  width = op.Width();
543 
544  // Initialize the processor grid if necessary
545  if (!gridInitialized_)
546  {
547  this->SetupGrid();
548  }
549 }
550 
551 } // mfem namespace
552 
553 #endif // MFEM_USE_MPI
554 #endif // MFEM_USE_SUPERLU
int Size() const
Logical size of the array.
Definition: array.hpp:110
void * InternalData() const
Definition: superlu.hpp:69
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: superlu.cpp:445
void * ScalePermstructPtr_
Definition: superlu.hpp:140
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
void SetRowPermutation(superlu::RowPerm row_perm, Array< int > *perm=NULL)
Definition: superlu.cpp:324
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
void SetSymmetricPattern(bool sym)
Definition: superlu.cpp:395
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: superlu.cpp:528
void SetIterativeRefine(superlu::IterRefine iter_ref)
Definition: superlu.cpp:361
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
void SetColumnPermutation(superlu::ColPerm col_perm)
Definition: superlu.cpp:315
void SetNumLookAheads(int num_lookaheads)
Definition: superlu.cpp:379
void SetPrintStatistics(bool print_stat)
Definition: superlu.cpp:297
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:69
void mfem_error(const char *msg)
Definition: error.cpp:107
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
void * SOLVEstructPtr_
Definition: superlu.hpp:142
void SetLookAheadElimTree(bool etree)
Definition: superlu.cpp:386
unsigned int sqrti(const unsigned int &a)
Definition: superlu.cpp:27
void SetReplaceTinyPivot(bool rtp)
Definition: superlu.cpp:370
Vector data type.
Definition: vector.hpp:41
void SetTranspose(superlu::Trans trans)
Definition: superlu.cpp:352
const SuperLURowLocMatrix * APtr_
Definition: superlu.hpp:130
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:54
Abstract operator.
Definition: operator.hpp:21
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
SuperLUSolver(MPI_Comm comm)
Definition: superlu.cpp:177
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
void SetEquilibriate(bool equil)
Definition: superlu.cpp:306