Amesos Package Browser (Single Doxygen Collection)  Development
Amesos_Superludist.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Amesos: Direct Sparse Solver Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 // TO DO: use Stat structure to print statistics???
30 // allow users to specify usermap ???
31 
32 #include "Amesos_Superludist.h"
33 #include "Epetra_Map.h"
34 #include "Epetra_Import.h"
35 #include "Epetra_CrsMatrix.h"
36 #include "Epetra_Util.h"
37 // #include "CrsMatrixTranspose.h"
38 #include "superlu_ddefs.h"
39 #include "supermatrix.h"
40 // SuperLU defines Reduce to be a macro in util.h, this conflicts with Reduce() in Epetra_MultiVector.h
41 #undef Reduce
42 
44 public:
45  // Teuchos::RCP<klu_symbolic> Symbolic_ ;
46  // Teuchos::RCP<klu_numeric> Numeric_ ;
47  fact_t FactOption_;
48  // Here are the structures used by Superlu
49  SuperMatrix SuperluA_;
50  ScalePermstruct_t ScalePermstruct_;
51  LUstruct_t LUstruct_;
52  SOLVEstruct_t SOLVEstruct_;
54  gridinfo_t grid_;
56 #if SUPERLU_DIST_MAJOR_VERSION == 5
57  //Note we may add the need for minor or patch version as need
58  superlu_dist_options_t options_;
59 #else
60  superlu_options_t options_;
61 #endif
62 
63 } ;
64 
65 
66 // ======================================================================
67 int Superludist_NumProcRows( int NumProcs ) {
68 #ifdef TFLOP
69  // Else, parameter 6 of DTRSV CTNLU is incorrect
70  return 1;
71 #else
72  int i;
73  int NumProcRows ;
74  for ( i = 1; i*i <= NumProcs; i++ )
75  ;
76  bool done = false ;
77  for ( NumProcRows = i-1 ; done == false ; ) {
78  int NumCols = NumProcs / NumProcRows ;
79  if ( NumProcRows * NumCols == NumProcs )
80  done = true;
81  else
82  NumProcRows-- ;
83  }
84  return NumProcRows;
85 #endif
86 }
87 
88 // ======================================================================
89 int SetNPRowAndCol(const int MaxProcesses, int& nprow, int& npcol)
90 {
91  nprow = Superludist_NumProcRows(MaxProcesses);
92  if (nprow < 1 ) nprow = 1;
93  npcol = MaxProcesses / nprow;
94 
95  if( nprow <=0 || npcol <= 0 || MaxProcesses <=0 ) {
96  std::cerr << "Amesos_Superludist: wrong value for MaxProcs ("
97  << MaxProcesses << "), or nprow (" << nprow
98  << ") or npcol (" << npcol << ")" << std::endl;
99  AMESOS_CHK_ERR(-1);
100  }
101  return(0);
102 }
103 
104 //=============================================================================
106  PrivateSuperluData_( rcp( new Amesos_Superlu_Pimpl() ) ),
107  Problem_(&prob),
108  RowMatrixA_(0),
109  GridCreated_(0),
110  FactorizationDone_(0),
111  FactorizationOK_(false),
112  NumGlobalRows_(0),
113  nprow_(0),
114  npcol_(0),
115  PrintNonzeros_(false),
116  ColPerm_("NOT SET"),
117  RowPerm_("NOT SET"),
118  perm_c_(0),
119  perm_r_(0),
120  IterRefine_("NOT SET"),
121  ReplaceTinyPivot_(true),
122  MtxConvTime_(-1),
123  MtxRedistTime_(-1),
124  VecRedistTime_(-1),
125  NumFactTime_(-1),
126  SolveTime_(-1),
127  OverheadTime_(-1)
128 {
129  Redistribute_ = true;
130  AddZeroToDiag_ = false;
131  PrivateSuperluData_->FactOption_ = SamePattern_SameRowPerm ;
132  ReuseSymbolic_ = false ;
133 
134  MaxProcesses_ = - 1;
135  Equil_ = true;
136  ColPerm_ = "MMD_AT_PLUS_A";
137  perm_c_ = 0;
138  RowPerm_ = "LargeDiag";
139  perm_r_ = 0;
140  IterRefine_ = "DOUBLE";
141  ReplaceTinyPivot_ = true;
142 
143  PrintNonzeros_ = false;
144 
145  ComputeTrueResidual_ = false;
146  ComputeVectorNorms_ = false;
147  PrintStatus_ = false ;
148  PrintTiming_ = false ;
149 
150  Teuchos::ParameterList ParamList;
151  SetParameters(ParamList);
152 }
153 
154 //=============================================================================
156 {
157  if (PrintTiming_) PrintTiming();
158  if (PrintStatus_) PrintStatus();
159 
160  if ( FactorizationDone_ ) {
161  SUPERLU_FREE( PrivateSuperluData_->SuperluA_.Store );
162  ScalePermstructFree(&PrivateSuperluData_->ScalePermstruct_);
164  LUstructFree(&PrivateSuperluData_->LUstruct_);
165  if ( PrivateSuperluData_->options_.SolveInitialized ) {
167  }
168  }
169  if ( GridCreated_ ) {
170  superlu_gridexit(&PrivateSuperluData_->grid_);
171  }
172 }
173 
174 // ======================================================================
176 {
177  // retrive general parameters
178 
180 
182 
183  if (ParameterList.isParameter("Redistribute"))
184  Redistribute_ = ParameterList.get<bool>("Redistribute");
185 
186  // parameters for Superludist only
187 
188  if (ParameterList.isSublist("Superludist") )
189  {
190  const Teuchos::ParameterList& SuperludistParams =
191  ParameterList.sublist("Superludist") ;
192 
193  if( SuperludistParams.isParameter("ReuseSymbolic") )
194  ReuseSymbolic_ = SuperludistParams.get<bool>("ReuseSymbolic");
195  std::string FactOption = "NotSet";
196 
197  if( SuperludistParams.isParameter("Fact") )
198  FactOption = SuperludistParams.get<std::string>("Fact");
199 
200  if( FactOption == "SamePattern_SameRowPerm" ) PrivateSuperluData_->FactOption_ = SamePattern_SameRowPerm;
201  else if( FactOption == "SamePattern" ) PrivateSuperluData_->FactOption_ = SamePattern;
202  else if ( FactOption != "NotSet" )
203  AMESOS_CHK_ERR(-2); // input not valid
204 
205  if (SuperludistParams.isParameter("Equil"))
206  Equil_ = SuperludistParams.get<bool>("Equil");
207 
208  if (SuperludistParams.isParameter("ColPerm"))
209  ColPerm_ = SuperludistParams.get<std::string>("ColPerm");
210 
211  if (ColPerm_ == "MY_PERMC")
212  if( SuperludistParams.isParameter("perm_c"))
213  perm_c_ = SuperludistParams.get<int*>("perm_c");
214 
215  if (SuperludistParams.isParameter("RowPerm"))
216  RowPerm_ = SuperludistParams.get<std::string>("RowPerm");
217  if( RowPerm_ == "MY_PERMR" ) {
218  if (SuperludistParams.isParameter("perm_r"))
219  perm_r_ = SuperludistParams.get<int*>("perm_r");
220  }
221 
222  if (SuperludistParams.isParameter("IterRefine"))
223  IterRefine_ = SuperludistParams.get<std::string>("IterRefine");
224 
225  if (SuperludistParams.isParameter("ReplaceTinyPivot"))
226  ReplaceTinyPivot_ = SuperludistParams.get<bool>("ReplaceTinyPivot");
227 
228  if (SuperludistParams.isParameter("PrintNonzeros"))
229  PrintNonzeros_ = SuperludistParams.get<bool>("PrintNonzeros");
230  }
231 
232  return(0);
233 }
234 
235 // ======================================================================
236 // Tasks of this method:
237 // 1) To set the required number of processes;
238 // 2) To set nprow_ and npcol_
239 // 3) To create a linear distribution (map) with elements only in the
240 // active processes
241 // 4) To redistribute the matrix from the original to the linear map.
242 // ======================================================================
244 {
245  ResetTimer(0);
246 
248  AMESOS_CHK_ERR(-1); // something has changed
249 
250  int iam = Comm().MyPID();
251 
254 
255  int m_per_p = NumGlobalRows_ / MaxProcesses_ ;
256  int remainder = NumGlobalRows_ - ( m_per_p * MaxProcesses_ );
257  int MyFirstElement = iam * m_per_p + EPETRA_MIN( iam, remainder );
258  int MyFirstNonElement = (iam+1) * m_per_p + EPETRA_MIN( iam+1, remainder );
259  int NumExpectedElements = MyFirstNonElement - MyFirstElement ;
260 
261  if (iam >= MaxProcesses_)
262  NumExpectedElements = 0;
263 
265  AMESOS_CHK_ERR(-1);
266 
267  const Epetra_Map &OriginalMap = RowMatrixA_->RowMatrixRowMap();
268 
269  UniformMap_ = rcp(new Epetra_Map(NumGlobalRows_, NumExpectedElements,
270  0, Comm()));
272  UniformMatrix_ = rcp(&CrsUniformMatrix(), false);
273  Importer_ = rcp(new Epetra_Import(*UniformMap_, OriginalMap));
274 
276  // int TracebackMode = Comm().GetTracebackMode();
277  // UniformMatrix_->SetTracebackMode(0);
278  if (AddZeroToDiag_ ) {
279  //
280  // Add 0.0 to each diagonal entry to avoid empty diagonal entries;
281  //
282  int OriginalTracebackMode = CrsUniformMatrix_->GetTracebackMode() ;
283  CrsUniformMatrix_->SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) ) ;
284  double zero = 0.0;
285  for ( int i = 0 ; i < UniformMap_->NumGlobalElements(); i++ )
286  if ( CrsUniformMatrix_->LRID(i) >= 0 )
287  CrsUniformMatrix_->InsertGlobalValues( i, 1, &zero, &i ) ;
288  CrsUniformMatrix_->SetTracebackMode( OriginalTracebackMode ) ;
289  }
290  // UniformMatrix_->SetTracebackMode(TracebackMode);
291 
292 
294 
295  MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
296 
297 
298 
299  return(0);
300 }
301 
302 
303 // ======================================================================
305 {
306  ResetTimer(1); // for "overhead"
307 
308  // FIXME????
309  // For now, if you change the shape of a matrix, you need to
310  // create a new Amesos instance.
311  //
312  //
314  AMESOS_CHK_ERR(-5);
315 
317 
318  if (Comm().NumProc() == 1)
319  Redistribute_ = false;
320 
321  // Set the matrix and grid shapes. Time is tracked within
322  // the RedistributeA() function
323 
324  if (Redistribute_)
325  RedistributeA() ;
326  else
327  {
328  if (Comm().NumProc() == 1)
329  {
330  nprow_ = 1;
331  npcol_ = 1;
332  }
333  else
334  {
336  AMESOS_CHK_ERR(-2);
337  SetNPRowAndCol(Comm().NumProc(), nprow_, npcol_);
338  }
339 
340  UniformMatrix_ = rcp(RowMatrixA_, false);
341  }
342 
343  // Extract Ai_, Ap_ and Aval_ from UniformMatrix_
344 
345  ResetTimer(0);
346 
347  int MyActualFirstElement = UniformMatrix().RowMatrixRowMap().MinMyGID() ;
348  int NumMyElements = UniformMatrix().NumMyRows() ;
349  int nnz_loc = UniformMatrix().NumMyNonzeros() ;
350  Ap_.resize( NumMyElements+1 );
351  Ai_.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ;
352  Aval_.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ;
353 
354  int NzThisRow ;
355  int Ai_index = 0 ;
356  int MyRow;
357  //int num_my_cols = UniformMatrix().NumMyCols() ;
358  double *RowValues;
359  int *ColIndices;
360  int MaxNumEntries_ = UniformMatrix().MaxNumEntries();
361  std::vector<double> RowValuesV_(MaxNumEntries_);
362  std::vector<int> ColIndicesV_(MaxNumEntries_);
363 
365 
366  const Epetra_CrsMatrix *SuperluCrs = dynamic_cast<const Epetra_CrsMatrix *>(&UniformMatrix());
367 
368  int ierr;
369 
370  for (MyRow = 0; MyRow < NumMyElements ; MyRow++)
371  {
372  if (SuperluCrs != 0)
373  {
374  ierr = SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues,
375  ColIndices);
376 
377  }
378  else {
379  ierr = UniformMatrix().ExtractMyRowCopy(MyRow, MaxNumEntries_, NzThisRow,
380  &RowValuesV_[0], &ColIndicesV_[0]);
381  RowValues = &RowValuesV_[0];
382  ColIndices = &ColIndicesV_[0];
383  }
384 
385  AMESOS_CHK_ERR(ierr);
386 
387  // MS // Added on 15-Mar-05
388  if (AddToDiag_ != 0.0 || AddZeroToDiag_) {
389  for (int i = 0 ; i < NzThisRow ; ++i) {
390  if (ColIndices[i] == MyRow) {
391  RowValues[i] += AddToDiag_;
392  break;
393  }
394  }
395  }
396 
397  Ap_[MyRow] = Ai_index ;
398  for ( int j = 0; j < NzThisRow; j++ ) {
399  Ai_[Ai_index] = Global_Columns_[ColIndices[j]] ;
400  Aval_[Ai_index] = RowValues[j] ;
401  Ai_index++;
402  }
403  }
404  assert( NumMyElements == MyRow );
405  Ap_[ NumMyElements ] = Ai_index ;
406 
407  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
408 
409  //
410  // Setup Superlu's grid
411  //
412  const Epetra_MpiComm & comm1 = dynamic_cast<const Epetra_MpiComm &> (Comm());
413 
414  if ( ! GridCreated_ ) {
415  // NOTE: nprow_ and npcol_ cannot be changed by the user
416  GridCreated_ = true;
417  superlu_gridinit(comm1.Comm(), nprow_, npcol_, &PrivateSuperluData_->grid_);
418  }
419 
420  if ( FactorizationDone_ ) {
421  SUPERLU_FREE( PrivateSuperluData_->SuperluA_.Store );
422  ScalePermstructFree(&PrivateSuperluData_->ScalePermstruct_);
424  LUstructFree(&PrivateSuperluData_->LUstruct_);
425  if ( PrivateSuperluData_->options_.SolveInitialized ) {
427  }
428  }
429 
430  MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
431  ResetTimer(0);
432 
433  //
434  // Only those processes in the grid participate from here on
435  //
436  if (Comm().MyPID() < nprow_ * npcol_) {
437  //
438  // Set up Superlu's data structures
439  //
440  set_default_options_dist(&PrivateSuperluData_->options_);
441 
442  dCreate_CompRowLoc_Matrix_dist( &PrivateSuperluData_->SuperluA_, NumGlobalRows_, NumGlobalRows_,
443  nnz_loc, NumMyElements, MyActualFirstElement,
444  &Aval_[0], &Ai_[0], &Ap_[0],
445  SLU_NR_loc, SLU_D, SLU_GE );
446 
447  FactorizationDone_ = true; // i.e. clean up Superlu data structures in the destructor
448 
450 #ifdef HAVE_SUPERLUDIST_LUSTRUCTINIT_2ARG
452 #else
454 #endif
455 
456  // stick options from ParameterList to options_ structure
457  // Here we follow the same order of the SuperLU_dist 2.0 manual (pag 55/56)
458 
459  assert( PrivateSuperluData_->options_.Fact == DOFACT );
460  PrivateSuperluData_->options_.Fact = DOFACT ;
461 
462  if( Equil_ ) PrivateSuperluData_->options_.Equil = (yes_no_t)YES;
463  else PrivateSuperluData_->options_.Equil = NO;
464 
465  if( ColPerm_ == "NATURAL" ) PrivateSuperluData_->options_.ColPerm = NATURAL;
466  else if( ColPerm_ == "MMD_AT_PLUS_A" ) PrivateSuperluData_->options_.ColPerm = MMD_AT_PLUS_A;
467  else if( ColPerm_ == "MMD_ATA" ) PrivateSuperluData_->options_.ColPerm = MMD_ATA;
468  // else if( ColPerm_ == "COLAMD" ) PrivateSuperluData_->options_.ColPerm = COLAMD; // COLAMD no longer supported in Superludist, as of July 2005
469  else if( ColPerm_ == "MY_PERMC" ) {
470  PrivateSuperluData_->options_.ColPerm = MY_PERMC;
472  }
473 
474  if( RowPerm_ == "NATURAL" ) PrivateSuperluData_->options_.RowPerm = (rowperm_t)NATURAL;
475  if( RowPerm_ == "LargeDiag" ) PrivateSuperluData_->options_.RowPerm = LargeDiag;
476  else if( ColPerm_ == "MY_PERMR" ) {
477  PrivateSuperluData_->options_.RowPerm = MY_PERMR;
479  }
480 
481  if( ReplaceTinyPivot_ ) PrivateSuperluData_->options_.ReplaceTinyPivot = (yes_no_t)YES;
482  else PrivateSuperluData_->options_.ReplaceTinyPivot = (yes_no_t)NO;
483 
484  if( IterRefine_ == "NO" ) PrivateSuperluData_->options_.IterRefine = (IterRefine_t)NO;
485  else if( IterRefine_ == "DOUBLE" ) {
486  PrivateSuperluData_->options_.IterRefine =
487 #ifdef HAVE_SUPERLUDIST_ENUM_NAMESPACE
488  SLU_DOUBLE
489 #else
490  DOUBLE
491 #endif
492  ;
493  }
494  else if( IterRefine_ == "EXTRA" ) {
495  PrivateSuperluData_->options_.IterRefine =
496 #ifdef HAVE_SUPERLUDIST_ENUM_NAMESPACE
497  SLU_EXTRA
498 #else
499  EXTRA
500 #endif
501  ;
502  }
503 
504  // Without the following two lines, SuperLU_DIST cannot be made
505  // quiet.
506 
507  if (PrintNonzeros_) PrivateSuperluData_->options_.PrintStat = (yes_no_t)YES;
508  else PrivateSuperluData_->options_.PrintStat = (yes_no_t)NO;
509 
510  SuperLUStat_t stat;
511  PStatInit(&stat); /* Initialize the statistics variables. */
512 
513  //
514  // Factor A using Superludsit (via a call to pdgssvx)
515  //
516  int info ;
517  double berr ; // Should be untouched
518  double xValues; // Should be untouched
519  int nrhs = 0 ; // Prevents forward and back solves
520  int ldx = NumGlobalRows_; // Should be untouched
521 
524 
525  if ( PrivateSuperluData_->options_.SolveInitialized ) {
527  }
528  AMESOS_CHK_ERR(info);
529 
530  PStatFree(&stat);
531  }
532 
533  NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
534 
535  return 0;
536 }
537 
538 // ======================================================================
539 // Refactor - Refactor the matrix
540 //
541 // Preconditions:
542 // The non-zero pattern of the matrix must not have changed since the
543 // previous call to Factor(). Refactor ensures that each process owns
544 // the same number of columns that it did on the previous call to Factor()
545 // and returns -4 if a discrepancy is found. However, that check does not
546 // guarantee that no change was made to the non-zero structure of the matrix.
547 // No call to SetParameters should be made between the call to Factor()
548 // and the call to Refactor(). If the user does not call SetParameters,
549 // as they need never do, they are safe on this.
550 //
551 // Postconditions:
552 // The matrix specified by Problem_->Operator() will have been redistributed,
553 // converted to the form needed by Superludist and factored.
554 // Ai_, Aval_
555 // SuperluA_
556 // SuperLU internal data structures reflecting the LU factorization
557 // ScalePermstruct_
558 // LUstructInit_
559 //
560 // Performance notes:
561 // Refactor does not allocate or de-allocate memory.
562 //
563 // Return codes:
564 // -4 if we detect a change to the non-zero structure of the matrix.
565 //
567 {
568  ResetTimer(0);
569  ResetTimer(1);
570 
571  //
572  // Update Ai_ and Aval_ (while double checking Ap_)
573  //
574  if (Redistribute_)
575  if(CrsUniformMatrix().Import(*RowMatrixA_, Importer(), Insert))
576  AMESOS_CHK_ERR(-4);
577 
578  MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_, 0);
579  ResetTimer(0);
580 
581  const Epetra_CrsMatrix *SuperluCrs = dynamic_cast<const Epetra_CrsMatrix *>(&UniformMatrix());
582 
583  double *RowValues;
584  int *ColIndices;
585  int MaxNumEntries_ = UniformMatrix().MaxNumEntries();
586  int NumMyElements = UniformMatrix().NumMyRows() ;
587  std::vector<int> ColIndicesV_(MaxNumEntries_);
588  std::vector<double> RowValuesV_(MaxNumEntries_);
589 
590  int NzThisRow ;
591  int Ai_index = 0 ;
592  int MyRow;
593  int ierr;
594 
595  for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
596  if ( SuperluCrs != 0 ) {
597  ierr = SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues,
598  ColIndices);
599  }
600  else {
601  ierr = UniformMatrix().ExtractMyRowCopy( MyRow, MaxNumEntries_,
602  NzThisRow, &RowValuesV_[0],
603  &ColIndicesV_[0]);
604  RowValues = &RowValuesV_[0];
605  ColIndices = &ColIndicesV_[0];
606  }
607 
608  AMESOS_CHK_ERR(ierr);
609 
610  if ( Ap_[MyRow] != Ai_index ) AMESOS_CHK_ERR(-4);
611  for ( int j = 0; j < NzThisRow; j++ ) {
612  // pdgssvx alters Ai_, so we have to set it again.
613  Ai_[Ai_index] = Global_Columns_[ColIndices[j]];
614  Aval_[Ai_index] = RowValues[j] ;
615  Ai_index++;
616  }
617  }
618  if( Ap_[ NumMyElements ] != Ai_index ) AMESOS_CHK_ERR(-4);
619 
620  MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
621  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
622  ResetTimer(0);
623 
624  if (Comm().MyPID() < nprow_ * npcol_) {
625 
626 
627  // If we reuse the same options, the code fails on multiprocess runs
628  set_default_options_dist(&PrivateSuperluData_->options_);
629 
630  if (PrintNonzeros_) PrivateSuperluData_->options_.PrintStat = (yes_no_t)YES;
631  else PrivateSuperluData_->options_.PrintStat = (yes_no_t)NO;
632 
633 
635  SuperLUStat_t stat;
636  PStatInit(&stat); /* Initialize the statistics variables. */
637  int info ;
638  double berr ; // Should be untouched
639  double xValues; // Should be untouched
640  int nrhs = 0 ; // Prevents forward and back solves
641  int ldx = NumGlobalRows_; // Should be untouched
644  PStatFree(&stat);
645  AMESOS_CHK_ERR( info ) ;
646  }
647 
648  NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
649 
650  return 0;
651 }
652 
653 // ======================================================================
655 {
656  if (GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
657  GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints())
658  return(false);
659  else
660  return(true);
661 }
662 
663 // ======================================================================
665 {
666  FactorizationOK_ = false ;
667 
668  return(0);
669 }
670 
671 // ======================================================================
673 {
674  RowMatrixA_ = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
675  if (RowMatrixA_ == 0)
676  AMESOS_CHK_ERR(-1); // Linear problem does not contain Epetra_RowMatrix
677 
678  // reset factorization
679  // FactorizationOK_ = false;
680 
681  if (!MatrixShapeOK())
682  AMESOS_CHK_ERR(-1); // matrix not square
683 
684  CreateTimer(Comm(), 2);
685 
687  ReFactor();
688  else
689  Factor();
690 
691  FactorizationOK_ = true;
692 
693  NumNumericFact_++;
694 
695  return(0);
696 }
697 
698 // ======================================================================
699 // We proceed as follows:
700 // - perform numeric factorization if not done;
701 // - extract solution and right-hand side;
702 // - redistribute them if necessary by creating additional
703 // working vectors, or use vecX and vecB otherwise;
704 // - call SuperLU_DIST's solve routine;
705 // - re-ship data if necessary.
706 // ======================================================================
708 {
709  if (!FactorizationOK_)
711 
712  ResetTimer(1);
713 
714  Epetra_MultiVector* vecX = Problem_->GetLHS();
715  Epetra_MultiVector* vecB = Problem_->GetRHS();
716 
717  if (vecX == 0 || vecB == 0)
718  AMESOS_CHK_ERR(-1);
719 
720  int nrhs = vecX->NumVectors() ;
721  if (vecB->NumVectors() != nrhs)
722  AMESOS_CHK_ERR(-1);
723 
724  double *values;
725  int ldx;
726 
727  RCP<Epetra_MultiVector> vec_uni;
728 
729  if (Redistribute_)
730  {
731  vec_uni = Teuchos::rcp(new Epetra_MultiVector(*UniformMap_, nrhs));
732  ResetTimer(0);
733  vec_uni->Import(*vecB, Importer(), Insert);
734  VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
735  }
736  else
737  {
738  vecX->Update(1.0, *vecB, 0.0);
739  vec_uni = Teuchos::rcp(vecX, false);
740  }
741 
742  //int NumMyElements = vec_uni->MyLength();
743  AMESOS_CHK_ERR(vec_uni->ExtractView(&values, &ldx));
744 
745  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
746 
747  ResetTimer(0);
748 
749  /* Bail out if I do not belong in the grid. */
750  if (Comm().MyPID() < nprow_ * npcol_)
751  {
752  int info ;
753  std::vector<double>berr(nrhs);
754  SuperLUStat_t stat;
755  PStatInit(&stat); /* Initialize the statistics variables. */
756 
758  AMESOS_CHK_ERR(-1); // internal error
759  PrivateSuperluData_->options_.Fact = FACTORED ;
760 
761  //bool BlockSolve = true ;
764  &stat, &info);
765  AMESOS_CHK_ERR(info);
766 
767  PStatFree(&stat);
768  }
769 
770  SolveTime_ = AddTime("Total solve time", SolveTime_, 0);
771 
772  ResetTimer(1);
773 
774  if (Redistribute_)
775  {
776  ResetTimer(0);
777  vecX->Export(*vec_uni, Importer(), Insert);
778  VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
779  }
780 
782  ComputeTrueResidual(*RowMatrixA_, *vecX, *vecB, UseTranspose(),
783  "Amesos_Superludist");
784 
786  ComputeVectorNorms(*vecX, *vecB, "Amesos_Superludist");
787 
788  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
789  NumSolve_++;
790 
791  return(0);
792 }
793 
794 // ======================================================================
796 {
797  if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
798  return;
799 
800  std::string p = "Amesos_Superludist : ";
801  int NNZ = RowMatrixA_->NumGlobalNonzeros();
802 
803  PrintLine();
804 
805  std::cout << p << "Matrix has " << NumGlobalRows_ << " rows"
806  << " and " << NNZ << " nonzeros" << std::endl;
807  std::cout << p << "Nonzero elements per row = "
808  << 1.0 * NNZ / NumGlobalRows_ << std::endl;
809  std::cout << p << "Percentage of nonzero elements = "
810  << 100.0 * NNZ /pow(NumGlobalRows_, 2.0) << std::endl;
811  std::cout << p << "Use transpose = " << UseTranspose() << std::endl;
812  std::cout << p << "Redistribute = " << Redistribute_ << std::endl;
813  std::cout << p << "# available processes = " << Comm().NumProc() << std::endl;
814  std::cout << p << "# processes used in computation = " << nprow_ * npcol_
815  << " ( = " << nprow_ << "x" << npcol_ << ")" << std::endl;
816  std::cout << p << "Equil = " << Equil_ << std::endl;
817  std::cout << p << "ColPerm = " << ColPerm_ << std::endl;
818  std::cout << p << "RowPerm = " << RowPerm_ << std::endl;
819  std::cout << p << "IterRefine = " << IterRefine_ << std::endl;
820  std::cout << p << "ReplaceTinyPivot = " << ReplaceTinyPivot_ << std::endl;
821  std::cout << p << "AddZeroToDiag = " << AddZeroToDiag_ << std::endl;
822  std::cout << p << "Redistribute = " << Redistribute_ << std::endl;
823 
824  PrintLine();
825 
826  return;
827 }
828 
829 // ======================================================================
831 {
832  if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
833  return;
834 
835  double ConTime = GetTime(MtxConvTime_);
836  double MatTime = GetTime(MtxRedistTime_);
837  double VecTime = GetTime(VecRedistTime_);
838  double NumTime = GetTime(NumFactTime_);
839  double SolTime = GetTime(SolveTime_);
840  double OveTime = GetTime(OverheadTime_);
841 
842  if (NumNumericFact_)
843  NumTime /= NumNumericFact_;
844 
845  if (NumSolve_)
846  SolTime /= NumSolve_;
847 
848  std::string p = "Amesos_Superludist : ";
849  PrintLine();
850 
851  std::cout << p << "Time to convert matrix to Superludist format = "
852  << ConTime << " (s)" << std::endl;
853  std::cout << p << "Time to redistribute matrix = "
854  << MatTime << " (s)" << std::endl;
855  std::cout << p << "Time to redistribute vectors = "
856  << VecTime << " (s)" << std::endl;
857  std::cout << p << "Number of symbolic factorizations = "
858  << NumSymbolicFact_ << std::endl;
859  std::cout << p << "Time for sym fact = 0.0 (s), avg = 0.0 (s)" << std::endl;
860  std::cout << p << "Number of numeric factorizations = "
861  << NumNumericFact_ << std::endl;
862  std::cout << p << "Time for num fact = "
863  << NumTime * NumNumericFact_ << " (s), avg = " << NumTime << " (s)" << std::endl;
864  std::cout << p << "Number of solve phases = "
865  << NumSolve_ << std::endl;
866  std::cout << p << "Time for solve = "
867  << SolTime * NumSolve_ << " (s), avg = " << SolTime << " (s)" << std::endl;
868 
869  double tt = NumTime * NumNumericFact_ + SolTime * NumSolve_;
870  if (tt != 0)
871  {
872  std::cout << p << "Total time spent in Amesos = " << tt << " (s) " << std::endl;
873  std::cout << p << "Total time spent in the Amesos interface = " << OveTime << " (s)" << std::endl;
874  std::cout << p << "(the above time does not include SuperLU_DIST time)" << std::endl;
875  std::cout << p << "Amesos interface time / total time = " << OveTime / tt << std::endl;
876  }
877 
878  PrintLine();
879 
880  return;
881 }
#define DOUBLE
void PrintTiming() const
Print various timig.
int NumSymbolicFact_
Number of symbolic factorization phases.
Definition: Amesos_Status.h:67
double GetTime(const std::string what) const
Gets the cumulative time using the string.
Definition: Amesos_Time.h:102
int MyGlobalElements(int *MyGlobalElementList) const
virtual const Epetra_Map & RowMatrixRowMap() const=0
int Solve()
Solves A X = B (or AT x = B)
bool ReuseSymbolic_
Allows FactOption to be used on subsequent calls to pdgssvx from NumericFactorization.
void PrintStatus() const
Print various information about the parameters used by Superludist.
MPI_Comm Comm() const
bool FactorizationOK_
true if NumericFactorization() has been successfully called.
int LRID(int GRID_in) const
std::vector< int > Ap_
T & get(ParameterList &l, const std::string &name)
bool UseTranspose() const
Always returns false. Superludist doesn&#39;t support transpose solve.
~Amesos_Superludist(void)
Amesos_Superludist Destructor.
bool Redistribute_
redistribute the input matrix prior to calling Superludist
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
static void SetTracebackMode(int TracebackModeValue)
bool isSublist(const std::string &name) const
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
ScalePermstruct_t ScalePermstruct_
bool AddZeroToDiag_
Adds zero to diagonal of redistributed matrix (some solvers choke on a matrix with a partly empty dia...
Insert
RCP< ParameterList > sublist(const RCP< ParameterList > &paramList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
void SetMaxProcesses(int &MaxProcesses, const Epetra_RowMatrix &A)
Definition: Amesos_Utils.h:77
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
Definition: Amesos_Time.h:64
#define EPETRA_MIN(x, y)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
int NumNumericFact_
Number of numeric factorization phases.
Definition: Amesos_Status.h:69
virtual int NumGlobalNonzeros() const=0
int GridCreated_
true if the SuperLU_DIST&#39;s grid has been created (and has to be free&#39;d)
virtual int MyPID() const=0
int FillComplete(bool OptimizeDataStorage=true)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int NumSolve_
Number of solves.
Definition: Amesos_Status.h:71
Epetra_CrsMatrix & CrsUniformMatrix()
Epetra_RowMatrix * RowMatrixA_
virtual int MaxNumEntries() const=0
std::vector< double > Aval_
const Epetra_Import & Importer() const
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
Definition: Amesos_Status.h:56
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
RCP< Epetra_Map > UniformMap_
#define AMESOS_CHK_ERR(a)
int NumVectors() const
void ComputeTrueResidual(const Epetra_RowMatrix &Matrix, const Epetra_MultiVector &X, const Epetra_MultiVector &B, const bool UseTranspose, const std::string prefix) const
Computes the true residual, B - Matrix * X, and prints the results.
Definition: Amesos_Utils.h:29
bool PrintTiming_
If true, prints timing information in the destructor.
Definition: Amesos_Status.h:52
virtual int NumMyRows() const=0
gridinfo_t grid_
SuperLU_DIST&#39;s grid information.
int * Global_Columns_
Contains the global ID of local columns.
bool PrintStatus_
If true, print additional information in the destructor.
Definition: Amesos_Status.h:54
Teuchos::RCP< Amesos_Superlu_Pimpl > PrivateSuperluData_
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
Definition: Amesos_Time.h:80
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Epetra_MultiVector * GetRHS() const
RCP< Epetra_RowMatrix > UniformMatrix_
virtual int NumProc() const=0
int NumGlobalRows_
Global dimension of the matrix.
int MinMyGID() const
int SetNPRowAndCol(const int MaxProcesses, int &nprow, int &npcol)
const Epetra_RowMatrix & UniformMatrix() const
static int GetTracebackMode()
Amesos_Superludist(const Epetra_LinearProblem &LinearProblem)
Amesos_Superludist Constructor.
int NumGlobalElements() const
Teuchos::RCP< Epetra_Import > Importer_
Copy
void ComputeVectorNorms(const Epetra_MultiVector &X, const Epetra_MultiVector &B, const std::string prefix) const
Computes the norms of X and B and print the results.
Definition: Amesos_Utils.h:52
bool isParameter(const std::string &name) const
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
Definition: Amesos_Time.h:74
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
Epetra_MultiVector * GetLHS() const
virtual const Epetra_Map & RowMatrixColMap() const=0
int ExtractView(double **A, int *MyLDA) const
int Superludist_NumProcRows(int NumProcs)
bool LinearMap() const
Add
Epetra_Operator * GetOperator() const
const Epetra_LinearProblem * Problem_
bool MatrixShapeOK() const
Returns true if SUPERLUDIST can handle this matrix shape.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
RCP< Epetra_CrsMatrix > CrsUniformMatrix_
virtual int NumGlobalRows() const=0
virtual int NumMyNonzeros() const=0
void PrintLine() const
Prints line on std::cout.
Definition: Amesos_Utils.h:71
#define EPETRA_MAX(x, y)
std::vector< int > Ai_
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
Definition: Amesos_Status.h:58
superlu_options_t options_
Vector of options.
double AddToDiag_
Add this value to the diagonal.