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