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