Amesos Package Browser (Single Doxygen Collection)  Development
Amesos_Superlu.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 #include "Amesos_Superlu.h"
30 #include "Epetra_Map.h"
31 #include "Epetra_Import.h"
32 #include "Epetra_RowMatrix.h"
33 #include "Epetra_CrsMatrix.h"
34 #include "Epetra_Vector.h"
35 #include "Epetra_Util.h"
36 #include "Epetra_Time.h"
37 #include "Epetra_Comm.h"
38 #include "Epetra_LinearProblem.h"
39 
40 namespace SLU {
41  extern "C" {
42 #ifdef AMESOS_SUPERLU_PRE_JULY2005
43 # include "dsp_defs.h"
44 #else
45 # include "slu_ddefs.h"
46 #endif
47  }
48 } // namespace SLU
49 
50 class SLUData {
51 public:
52  SLU::SuperMatrix A, B, X, L, U;
53 #ifdef USE_DGSTRF
54  SLU::SuperMatrix AC;
55 #endif
56  SLU::superlu_options_t SLU_options;
57  SLU::mem_usage_t mem_usage;
58  SLU::fact_t refactor_option ; // SamePattern or SamePattern_SameRowPerm
59 
60  SLUData() {
61  A.Store = B.Store = X.Store = L.Store = U.Store = NULL;
62 #ifdef USE_DGSTRF
63  AC.Store = NULL;
64 #endif
65  SLU::set_default_options(&SLU_options);
66  refactor_option = SLU::DOFACT;
67  }
68 };
69 
70 using namespace Teuchos;
71 
72 //=============================================================================
74  DummyArray(NULL),
75  NumGlobalRows_(-1),
76  NumGlobalNonzeros_(-1),
77  UseTranspose_(false),
78  FactorizationOK_(false),
79  FactorizationDone_(false),
80  iam_(0),
81  MtxConvTime_(-1),
82  MtxRedistTime_(-1),
83  VecRedistTime_(-1),
84  NumFactTime_(-1),
85  SolveTime_(-1),
86  OverheadTime_(-1),
87  SerialMap_(Teuchos::null),
88  SerialCrsMatrixA_(Teuchos::null),
89  ImportToSerial_(Teuchos::null),
90  SerialMatrix_(0),
91  RowMatrixA_(0)
92 {
93  data_ = new SLUData();
94  ferr_.resize(1);
95  berr_.resize(1);
96 
97  Problem_ = &prob ;
98 
99  // I have to skip timing on this, because the Comm object is not done yet
100  dCreate_Dense_Matrix( &(data_->X),
101  0,
102  0,
103  DummyArray,
104  0,
105  SLU::SLU_DN, SLU::SLU_D, SLU::SLU_GE);
106 
107  dCreate_Dense_Matrix( &(data_->B),
108  0,
109  0,
110  DummyArray,
111  0,
112  SLU::SLU_DN, SLU::SLU_D, SLU::SLU_GE);
113 }
114 
115 //=============================================================================
117 {
118  if (PrintTiming_) PrintTiming();
119  if (PrintStatus_) PrintStatus();
120 
121  Destroy_SuperMatrix_Store(&data_->B);
122  Destroy_SuperMatrix_Store(&data_->X);
123 
124  if (iam_ == 0) {
125  if (FactorizationDone_) {
126  Destroy_SuperMatrix_Store(&data_->A);
127  Destroy_SuperNode_Matrix(&data_->L);
128  Destroy_CompCol_Matrix(&data_->U);
129  }
130  }
131 
132  delete data_;
133 }
134 
135 // ======================================================================
137 {
138  // retrive general parameters
139 
141 
143 
144  /* the list is empty now
145  if (ParameterList.isSublist("Superlu") ) {
146  Teuchos::ParameterList SuperluParams = ParameterList.sublist("Superlu") ;
147  }
148  */
149  return 0;
150 }
151 
152 // ======================================================================
154 {
155  if (GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints64() !=
156  GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints64())
157  return(false);
158 
159  return(true);
160 }
161 
162 // ======================================================================
164 {
165  ResetTimer(0);
166  ResetTimer(1);
167 
168  RowMatrixA_ = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
169  if (RowMatrixA_ == 0)
170  AMESOS_CHK_ERR(-1); // cannot cast to RowMatrix
171 
172  iam_ = Comm().MyPID() ;
173 
174  const Epetra_Map &OriginalMap = RowMatrixA_->RowMatrixRowMap() ;
175 
179  AMESOS_CHK_ERR(-1); // only square matrices
180  if (NumGlobalNonzeros_ <= 0)
181  AMESOS_CHK_ERR(-2); // empty matrix
182 
183  int NumMyElements_ = 0;
184  if (iam_ == 0) NumMyElements_ = NumGlobalRows_;
185 
186  // If the matrix is distributed, then brings it to processor zero.
187  // This also requires the construction of a serial map, and
188  // the exporter from distributed to serial.
189  // Otherwise, simply take the pointer of RowMatrixA_ and
190  // set it to SerialMatrix_.
191 
192  if (Comm().NumProc() == 1) // Bug #1411 - should recognize serial matrices even when NumProc() > 1
193  {
195  }
196  else
197  {
198 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES)
200  SerialMap_ = rcp(new Epetra_Map((int) NumGlobalRows_, NumMyElements_, 0, Comm()));
201  else
202 #endif
203 #if !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
205  SerialMap_ = rcp(new Epetra_Map((long long) NumGlobalRows_, NumMyElements_, 0, Comm()));
206  else
207 #endif
208  throw "Amesos_Superlu::ConvertToSerial: Global indices unknown.";
209 
210  ImportToSerial_ = rcp(new Epetra_Import(SerialMap(), OriginalMap));
211 
214 
215  // MS // Set zero element if not present, possibly add
216  // MS // something to the diagonal
217  //double AddToDiag = 0.0;
218  // FIXME?? bug #1371
219 #if 0
220  if (iam_ == 0)
221  {
222  int this_res ;
223  for (int i = 0 ; i < NumGlobalRows_; i++ ) {
224  // I am not sure what the intent is here,
225  // but InsertGlobalValues returns 1 meaning "out of room"
226  // and as a result, we sum AddToDiag_ into this diagonal element
227  // a second time.
228  if (this_res=SerialCrsMatrixA_->InsertGlobalValues(i, 1, &AddToDiag, &i)) {
229  std::cout << __FILE__ << "::" << __LINE__
230  << " this_res = " << this_res
231  << std::endl ;
232  SerialCrsMatrixA_->SumIntoGlobalValues(i, 1, &AddToDiag, &i);
233  }
234  }
235  }
236 #endif
237 
240  }
241 
242  MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_, 0);
243  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
244 
245  return(0);
246 }
247 
248 //
249 // See also pre and post conditions in Amesos_Superlu.h
250 // Preconditions:
251 // SerialMatrix_ points to the matrix to be factored and solved
252 // NumGlobalRows_ has been set to the dimension of the matrix
253 // NumGlobalNonzeros_ has been set to the number of non-zeros in the matrix
254 // (i.e. ConvertToSerial() has been called)
255 // FactorizationDone_ and FactorizationOK_ must be accurate
256 //
257 // Postconditions:
258 // data->A
259 // FactorizationDone_ = true
260 //
261 // Issues:
262 //
263 //
265 {
266  // Convert matrix to the form that Superlu expects (Ap, Ai, Aval)
267  // I suppose that the matrix has already been formed on processor 0
268  // as SerialMatrix_.
269 
270  ResetTimer(0);
271  if (iam_ == 0)
272  {
273  ResetTimer(1);
274 
280  {
281  AMESOS_CHK_ERR(-1); // something fishy here
282  }
283 
285  Ap_.resize(NumGlobalRows_ + 1, 0);
288 
289  int NzThisRow ;
290  int Ai_index = 0 ;
291  int MyRow;
292  int MaxNumEntries_ = SerialMatrix_->MaxNumEntries();
293 
294  for (MyRow = 0; MyRow < NumGlobalRows_ ; MyRow++ )
295  {
296  Ap_[MyRow] = Ai_index;
297  int ierr;
298  ierr = SerialMatrix_->ExtractMyRowCopy(MyRow, MaxNumEntries_, NzThisRow,
299  &Aval_[Ai_index], &Ai_[Ai_index]);
300  AMESOS_CHK_ERR(ierr);
301  Ai_index += NzThisRow;
302  }
303 
304  Ap_[NumGlobalRows_] = Ai_index;
305 
306  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
307 
308  if ( FactorizationDone_ ) {
309  Destroy_SuperMatrix_Store(&data_->A);
310  Destroy_SuperNode_Matrix(&data_->L);
311  Destroy_CompCol_Matrix(&data_->U);
312  }
313 
314  /* Create matrix A in the format expected by SuperLU. */
315  dCreate_CompCol_Matrix( &(data_->A), NumGlobalRows_, NumGlobalRows_,
317  &Ai_[0], &Ap_[0], SLU::SLU_NR, SLU::SLU_D, SLU::SLU_GE );
318  }
319 
320  FactorizationDone_ = true;
321  MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
322 
323  return 0;
324 }
325 
326 // ======================================================================
327 // MS // What is the difference between ReFactor() and Factor()?
328 // ======================================================================
330 {
331  // Convert matrix to the form that Superlu expects (Ap, Ai, Aval)
332  // I suppose that ConvertToSerial() has already been called,
333  // there I have SerialMatrix_ stored at it should.
334  // Only processor 0 does something useful here.
335 
336  if (iam_ == 0)
337  {
338  ResetTimer(1);
339 
345  {
346  AMESOS_CHK_ERR(-1);
347  }
348 
349  Ap_.resize(NumGlobalRows_+ 1, 0);
352 
353  int NzThisRow ;
354  int Ai_index = 0 ;
355  int MyRow;
356  double *RowValues;
357  int *ColIndices;
358  std::vector<int> ColIndicesV_;
359  std::vector<double> RowValuesV_;
360  int MaxNumEntries_ = SerialMatrix_->MaxNumEntries();
361 
362  Epetra_CrsMatrix *SuperluCrs = dynamic_cast<Epetra_CrsMatrix *>(SerialMatrix_);
363  if ( SuperluCrs == 0 ) {
364  ColIndicesV_.resize(MaxNumEntries_);
365  RowValuesV_.resize(MaxNumEntries_);
366  }
367  for ( MyRow = 0; MyRow < NumGlobalRows_ ; MyRow++ ) {
368  if ( SuperluCrs != 0 ) {
369  AMESOS_CHK_ERR(SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues, ColIndices));
370  }
371  else {
372  AMESOS_CHK_ERR(SerialMatrix_->ExtractMyRowCopy(MyRow, MaxNumEntries_,NzThisRow, &RowValuesV_[0],
373  &ColIndicesV_[0]));
374  RowValues = &RowValuesV_[0];
375  ColIndices = &ColIndicesV_[0];
376  }
377 
378  if (Ap_[MyRow] != Ai_index)
379  AMESOS_CHK_ERR(-4);
380 
381  for (int j = 0; j < NzThisRow; j++) {
382  assert(Ai_[Ai_index] == ColIndices[j]) ; // FIXME this may not work.
383  Aval_[Ai_index] = RowValues[j] ;
384  Ai_index++;
385  }
386  }
387  assert( NumGlobalRows_ == MyRow );
388  Ap_[ NumGlobalRows_ ] = Ai_index ;
389 
390  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
391 
392  assert ( FactorizationDone_ ) ;
393  Destroy_SuperMatrix_Store(&data_->A);
394  Destroy_SuperNode_Matrix(&data_->L);
395  Destroy_CompCol_Matrix(&data_->U);
396  /* Create matrix A in the format expected by SuperLU. */
397  dCreate_CompCol_Matrix( &(data_->A), NumGlobalRows_, NumGlobalRows_,
399  &Ai_[0], &Ap_[0], SLU::SLU_NR, SLU::SLU_D, SLU::SLU_GE );
400  }
401 
402  return 0;
403 }
404 
405 // ======================================================================
407 {
408  // nothing happens here, all it is done in NumericFactorization()
409  FactorizationOK_ = false;
410  return 0;
411 }
412 
413 // ======================================================================
415 {
416  CreateTimer(Comm(), 2);
417 
418  ResetTimer(0);
419 
420  ConvertToSerial();
421 
422  perm_r_.resize(NumGlobalRows_);
423  perm_c_.resize(NumGlobalRows_);
424  etree_.resize(NumGlobalRows_);
425  R_.resize(NumGlobalRows_);
426  C_.resize(NumGlobalRows_);
427 
428  SLU::superlu_options_t& SLUopt = data_->SLU_options ;
429 
430  set_default_options( &SLUopt ) ;
431  if (FactorizationOK_) {
433  SLUopt.Fact = data_->refactor_option ;
434  } else {
436  FactorizationOK_ = true;
437  SLUopt.Fact = SLU::DOFACT;
438  }
439 
440  int Ierr[1];
441  if ( iam_ == 0 ) {
442  double rpg, rcond;
443  equed_ = 'N';
444 
445 #if 0
446  if ( ! UseTranspose() ) // FIXME - I doubt we need this here.
447  assert( SLUopt.Trans == NOTRANS ) ;
448  else
449  SLUopt.Trans = TRANS ;
450 
451 
452  // SLUopt.ColPerm = COLAMD ;
453 
454  std::cout << " SLUopt.ColPerm = " << SLUopt.ColPerm << std::endl ;
455  std::cout << " SLUopt.Equil = " << SLUopt.Equil << std::endl ;
456  std::cout << " SLUopt.Fact = " << SLUopt.Fact << std::endl ;
457  std::cout << " SLUopt.IterRefine = " << SLUopt.IterRefine << std::endl ;
458  std::cout << " data_->A.Stype = " << data_->A.Stype
459  << " SLU_NC = " << SLU_NC
460  << " SLU_NR = " << SLU_NR
461  << std::endl ;
462  std::cout << " SLUopt.ColPerm = " << SLUopt.ColPerm << std::endl ;
463 #endif
464 
465  data_->B.nrow = NumGlobalRows_;
466  data_->B.ncol = 0;
467  SLU::DNformat* Bstore = (SLU::DNformat *) (data_->B.Store) ;
468  Bstore->lda = NumGlobalRows_;
469  Bstore->nzval = DummyArray;
470  data_->X.nrow = NumGlobalRows_;
471  data_->X.ncol = 0;
472  SLU::DNformat* Xstore = (SLU::DNformat *) (data_->X.Store) ;
473  Xstore->lda = NumGlobalRows_;
474  Xstore->nzval = DummyArray;
475 
476  SLU::SuperLUStat_t SLU_stat ;
477  SLU::StatInit( &SLU_stat ) ;
478  assert( SLUopt.Fact == SLU::DOFACT);
479  dgssvx( &(SLUopt), &(data_->A),
480  &perm_c_[0], &perm_r_[0], &etree_[0], &equed_, &R_[0],
481  &C_[0], &(data_->L), &(data_->U), NULL, 0,
482  &(data_->B), &(data_->X), &rpg, &rcond, &ferr_[0],
483  &berr_[0], &(data_->mem_usage), &SLU_stat,
484  &Ierr[0] );
485  SLU::StatFree( &SLU_stat ) ;
486  }
487 
488  NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
489 
490  FactorizationDone_ = true;
491 
492  ++NumNumericFact_;
493 
494  return(0);
495 }
496 
497 // ======================================================================
499 {
500  if (!FactorizationDone_) {
501  FactorizationOK_ = false;
503  }
504 
505  ResetTimer(1); // for "overhead'
506 
509  int Ierr;
510 
511  if (vecX == 0 || vecB == 0)
512  AMESOS_CHK_ERR(-1); // vectors not set
513 
514  int nrhs = vecX->NumVectors();
515  if (nrhs != vecB->NumVectors())
516  AMESOS_CHK_ERR(-2); // vectors not compatible
517 
518  ferr_.resize(nrhs);
519  berr_.resize(nrhs);
520 
521  Epetra_MultiVector* SerialB = 0;
522  Epetra_MultiVector* SerialX = 0;
523 
524  double *SerialXvalues ;
525  double *SerialBvalues ;
526 
527  if (Comm().NumProc() == 1)
528  {
529  SerialB = vecB;
530  SerialX = vecX;
531  }
532  else
533  {
534  ResetTimer(0);
535 
536  SerialX = new Epetra_MultiVector(SerialMap(), nrhs);
537  SerialB = new Epetra_MultiVector(SerialMap(), nrhs);
538 
539  SerialB->Import(*vecB, ImportToSerial(), Insert);
540  // SerialB = SerialB;
541  // SerialX = SerialX;
542 
543  VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
544  }
545 
546  ResetTimer(0);
547 
548  // Call SUPERLU's dgssvx to perform the solve
549  // At this point I have, on processor 0, the solution vector
550  // in SerialX and the right-hand side in SerialB
551 
552  if (iam_ == 0)
553  {
554  int SerialXlda;
555  int SerialBlda;
556 
557  int ierr;
558  ierr = SerialX->ExtractView(&SerialXvalues, &SerialXlda);
559  assert (ierr == 0);
560  assert (SerialXlda == NumGlobalRows_ ) ;
561 
562  ierr = SerialB->ExtractView(&SerialBvalues, &SerialBlda);
563  assert (ierr == 0);
564  assert (SerialBlda == NumGlobalRows_ ) ;
565 
566  SLU::SuperMatrix& dataX = (data_->X) ;
567  dataX.nrow = NumGlobalRows_;
568  dataX.ncol = nrhs ;
569  data_->X.nrow = NumGlobalRows_;
570  SLU::DNformat* Xstore = (SLU::DNformat *) (data_->X.Store) ;
571  Xstore->lda = SerialXlda;
572  Xstore->nzval = &SerialXvalues[0];
573 
574  SLU::SuperMatrix& dataB = (data_->B) ;
575  dataB.nrow = NumGlobalRows_;
576  dataB.ncol = nrhs ;
577  data_->B.nrow = NumGlobalRows_;
578  SLU::DNformat* Bstore = (SLU::DNformat *) (data_->B.Store) ;
579  Bstore->lda = SerialBlda;
580  Bstore->nzval = &SerialBvalues[0];
581 
582  double rpg, rcond;
583 #if 0
584  char fact, trans, refact;
585  fact = 'F';
586  refact = 'N';
587  trans = 'N' ; // This will allow us to try trans and not trans - see if either works.
588  // equed = 'N';
589 #endif
590 #if 0
591  dgssvx( &fact, &trans, &refact, &(data_->A), &(data_->iparam), &perm_c_[0],
592  &perm_r_[0], &etree_[0], &equed_, &R_[0], &C_[0], &(data_->L), &(data_->U),
593  NULL, 0, &(data_->B), &(data_->X), &rpg, &rcond,
594  &ferr_[0], &berr_[0], &(data_->mem_usage), &Ierr[0] );
595 #endif
596 
597  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1); // NOTE: only timings on processor 0 will be meaningful
598 
599  SLU::SuperLUStat_t SLU_stat ;
600  SLU::StatInit( &SLU_stat ) ;// Copy the scheme used in dgssvx1.c
601  data_->SLU_options.Fact = SLU::FACTORED ;
602  SLU::superlu_options_t& SLUopt = data_->SLU_options ;
603  if (UseTranspose())
604  SLUopt.Trans = SLU::TRANS;
605  else
606  SLUopt.Trans = SLU::NOTRANS;
607 
608  //#ifdef USE_DGSTRF
609  // assert( equed_ == 'N' ) ;
610  dgssvx( &(SLUopt), &(data_->A),
611  &perm_c_[0], &perm_r_[0], &etree_[0], &equed_, &R_[0],
612  &C_[0], &(data_->L), &(data_->U), NULL, 0,
613  &(data_->B), &(data_->X), &rpg, &rcond, &ferr_[0],
614  &berr_[0], &(data_->mem_usage), &SLU_stat,
615  &Ierr);
616  // assert( equed_ == 'N' ) ;
617  StatFree( &SLU_stat ) ;
618  }
619 
620  SolveTime_ = AddTime("Total solve time", SolveTime_, 0);
621 
622  ResetTimer(1); // for "overhead'
623 
624  //
625  // Copy X back to the original vector
626  //
627  if (Comm().NumProc() != 1)
628  {
629  ResetTimer(0);
630 
631  vecX->Export(*SerialX, ImportToSerial(), Insert);
632  delete SerialB;
633  delete SerialX;
634 
635  VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
636  }
637 
639  ComputeTrueResidual(*(GetProblem()->GetMatrix()), *vecX, *vecB,
640  UseTranspose(), "Amesos_Superlu");
641 
643  ComputeVectorNorms(*vecX, *vecB, "Amesos_Superlu");
644 
645  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
646 
647  ++NumSolve_;
648 
649  // All processes should return the same error code
650  if (Comm().NumProc() != 1)
651  Comm().Broadcast(&Ierr, 1, 0);
652 
653  AMESOS_RETURN(Ierr);
654 }
655 
656 // ================================================ ====== ==== ==== == =
657 
659 {
660  if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
661  return;
662 
663  std::string p = "Amesos_Superlu : ";
664  PrintLine();
665 
666  long long n = GetProblem()->GetMatrix()->NumGlobalRows64();
667  long long nnz = GetProblem()->GetMatrix()->NumGlobalNonzeros64();
668 
669  std::cout << p << "Matrix has " << n << " rows"
670  << " and " << nnz << " nonzeros" << std::endl;
671  std::cout << p << "Nonzero elements per row = "
672  << 1.0 * nnz / n << std::endl;
673  std::cout << p << "Percentage of nonzero elements = "
674  << 100.0 * nnz /(pow(double(n),double(2.0))) << std::endl;
675  std::cout << p << "Use transpose = " << UseTranspose_ << std::endl;
676 
677  PrintLine();
678 
679  return;
680 }
681 
682 // ================================================ ====== ==== ==== == =
684 {
685  if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
686  return;
687 
688  double ConTime = GetTime(MtxConvTime_);
689  double MatTime = GetTime(MtxRedistTime_);
690  double VecTime = GetTime(VecRedistTime_);
691  double NumTime = GetTime(NumFactTime_);
692  double SolTime = GetTime(SolveTime_);
693  double OveTime = GetTime(OverheadTime_);
694 
695  if (NumNumericFact_)
696  NumTime /= NumNumericFact_;
697 
698  if (NumSolve_)
699  SolTime /= NumSolve_;
700 
701  std::string p = "Amesos_Superlu : ";
702  PrintLine();
703 
704  std::cout << p << "Time to convert matrix to SuperLU format = " << ConTime << " (s)" << std::endl;
705  std::cout << p << "Time to redistribute matrix = " << MatTime << " (s)" << std::endl;
706  std::cout << p << "Time to redistribute vectors = " << VecTime << " (s)" << std::endl;
707  std::cout << p << "Number of numeric factorizations = " << NumNumericFact_ << std::endl;
708  std::cout << p << "Time for num fact = "
709  << NumTime * NumNumericFact_ << " (s), avg = "
710  << NumTime << " (s)" << std::endl;
711  std::cout << p << "Number of solve phases = " << NumSolve_ << std::endl;
712  std::cout << p << "Time for solve = "
713  << SolTime * NumSolve_ << " (s), avg = " << SolTime << " (s)" << std::endl;
714  double tt = NumTime * NumNumericFact_ + SolTime * NumSolve_;
715  if (tt != 0)
716  {
717  std::cout << p << "Total time spent in Amesos = " << tt << " (s) " << std::endl;
718  std::cout << p << "Total time spent in the Amesos interface = " << OveTime << " (s)" << std::endl;
719  std::cout << p << "(the above time does not include SuperLU time)" << std::endl;
720  std::cout << p << "Amesos interface time / total time = " << OveTime / tt << std::endl;
721  }
722 
723 
724  PrintLine();
725 
726  return;
727 }
double GetTime(const std::string what) const
Gets the cumulative time using the string.
Definition: Amesos_Time.h:102
std::vector< int > etree_
virtual const Epetra_Map & RowMatrixRowMap() const=0
Teuchos::RCP< Epetra_Map > SerialMap_
Contains a map with all elements assigned to processor 0.
std::vector< double > Aval_
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
UF_long CHOLMOD() nnz(cholmod_sparse *A, cholmod_common *Common)
Epetra_RowMatrix * SerialMatrix_
For parallel runs, stores the matrix defined on SerialMap_.
SLU::SuperMatrix B
SLU::fact_t refactor_option
SLU::superlu_options_t SLU_options
int Solve()
Solves A X = B (or AT x = B)
void PrintTiming() const
Prints timing information.
bool MatrixShapeOK() const
Returns true if the solver can handle this matrix shape.
std::vector< int > perm_c_
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
std::vector< int > Ap_
stores the matrix in SuperLU format.
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
SLU::SuperMatrix X
std::vector< double > C_
SLU::SuperMatrix L
SLU::SuperMatrix A
Insert
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
Definition: Amesos_Time.h:64
int Factor()
Factors the matrix, no previous factorization available.
int MtxConvTime_
Quick access pointer to internal timing data.
bool UseTranspose_
If true, solve the linear system with the transpose of the matrix.
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
long long NumGlobalNonzeros_
Global number of nonzeros in the matrix.
int NumNumericFact_
Number of numeric factorization phases.
Definition: Amesos_Status.h:69
T * get() const
virtual int MyPID() const=0
SLU::mem_usage_t mem_usage
int ReFactor()
Re-factors the matrix.
int FillComplete(bool OptimizeDataStorage=true)
virtual int NumMyCols() const=0
SLU::factor_param_t iparam
Definition: Epetra_SLU.cpp:49
int NumSolve_
Number of solves.
Definition: Amesos_Status.h:71
#define NULL
virtual int MaxNumEntries() const=0
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Amesos_Superlu(const Epetra_LinearProblem &LinearProblem)
Amesos_Superlu Constructor.
const Epetra_Import & ImportToSerial() const
Returns a reference to the importer.
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
Definition: Amesos_Status.h:56
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
#define AMESOS_CHK_ERR(a)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int NumVectors() const
bool GlobalIndicesInt() const
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
SLU::SuperMatrix U
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
std::vector< double > ferr_
bool UseTranspose() const
Returns the current UseTranspose setting.
bool PrintStatus_
If true, print additional information in the destructor.
Definition: Amesos_Status.h:54
#define AMESOS_RETURN(amesos_err)
TRANS
virtual long long NumGlobalCols64() const=0
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
virtual int Broadcast(double *MyVals, int Count, int Root) const=0
void PrintStatus() const
Prints status information.
Epetra_RowMatrix * RowMatrixA_
Pointer to the linear system matrix.
virtual long long NumGlobalNonzeros64() const=0
double * DummyArray
stores the matrix in SuperLU format.
int ConvertToSerial()
Sets up the matrix on processor 0.
Epetra_MultiVector * GetRHS() const
std::vector< double > R_
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrixA_
Contains a matrix with all rows assigned to processor 0.
Teuchos::RCP< Epetra_Import > ImportToSerial_
Importer from distributed to SerialMap_.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
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
SLUData * data_
Main structure for SuperLU.
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
Definition: Amesos_Time.h:74
~Amesos_Superlu()
Amesos_Superlu Destructor.
std::vector< double > berr_
std::vector< int > perm_r_
bool FactorizationOK_
If true, the factorization has been successfully computed.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
Epetra_RowMatrix * GetMatrix() const
Epetra_MultiVector * GetLHS() const
int ExtractView(double **A, int *MyLDA) const
int iam_
Process number (i.e. Comm().MyPID()
const Epetra_Map & SerialMap() const
Returns a reference to the serial map.
bool GlobalIndicesLongLong() const
Add
Epetra_Operator * GetOperator() const
long long NumGlobalRows_
Global size of the matrix.
int n
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
void PrintLine() const
Prints line on std::cout.
Definition: Amesos_Utils.h:71
std::vector< int > Ai_
stores the matrix in SuperLU format.
#define EPETRA_MAX(x, y)
virtual long long NumGlobalRows64() const=0
const Epetra_LinearProblem * Problem_
Pointer to the user&#39;s defined linear problem.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
Definition: Amesos_Status.h:58