Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Amesos_Klu.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_Klu.h"
30#include "Epetra_Map.h"
31#include "Epetra_Import.h"
32#include "Epetra_CrsMatrix.h"
33#include "Epetra_Vector.h"
34#include "Epetra_Util.h"
35#include "Amesos_Support.h"
36
37extern "C" {
38
39#include "trilinos_klu_decl.h"
40
41}
42
43namespace {
44
45using Teuchos::RCP;
46
47template<class T, class DeleteFunctor>
48class DeallocFunctorDeleteWithCommon
49{
50public:
51 DeallocFunctorDeleteWithCommon(
52 const RCP<trilinos_klu_common> &common
53 ,DeleteFunctor deleteFunctor
54 )
55 : common_(common), deleteFunctor_(deleteFunctor)
56 {}
57 typedef T ptr_t;
58 void free( T* ptr ) {
59 if(ptr) deleteFunctor_(&ptr,&*common_);
60 }
61private:
62 Teuchos::RCP<trilinos_klu_common> common_;
63 DeleteFunctor deleteFunctor_;
64 DeallocFunctorDeleteWithCommon(); // Not defined and not to be called!
65};
66
67template<class T, class DeleteFunctor>
68DeallocFunctorDeleteWithCommon<T,DeleteFunctor>
69deallocFunctorDeleteWithCommon(
70 const RCP<trilinos_klu_common> &common
71 ,DeleteFunctor deleteFunctor
72 )
73{
74 return DeallocFunctorDeleteWithCommon<T,DeleteFunctor>(common,deleteFunctor);
75}
76
77
78} // namespace
79
81{
82public:
83 Teuchos::RCP<trilinos_klu_symbolic> Symbolic_;
84 Teuchos::RCP<trilinos_klu_numeric> Numeric_;
85 Teuchos::RCP<trilinos_klu_common> common_;
86};
87
88//=============================================================================
89Amesos_Klu::Amesos_Klu(const Epetra_LinearProblem &prob ) :
90 PrivateKluData_( rcp( new Amesos_Klu_Pimpl() ) ),
91 CrsMatrixA_(0),
92 TrustMe_(false),
93 UseTranspose_(false),
94 Problem_(&prob),
95 MtxRedistTime_(-1),
96 MtxConvTime_(-1),
97 VecRedistTime_(-1),
98 SymFactTime_(-1),
99 NumFactTime_(-1),
100 SolveTime_(-1),
101 OverheadTime_(-1)
102{
103 // MS // move declaration of Problem_ above because I need it
104 // MS // set up before calling Comm()
105 Teuchos::ParameterList ParamList ;
106 SetParameters( ParamList ) ;
107
108}
109
110//=============================================================================
112
113 // print out some information if required by the user
114 if( (verbose_ && PrintTiming_) || verbose_ == 2 ) PrintTiming();
115 if( (verbose_ && PrintStatus_) || verbose_ == 2 ) PrintStatus();
116}
117
118//=============================================================================
120{
121
122 if ( AddZeroToDiag_==0 && numentries_ != RowMatrixA_->NumGlobalNonzeros64()) {
123 std::cerr << " The number of non zero entries in the matrix has changed since the last call to SymbolicFactorization(). " ;
124 AMESOS_CHK_ERR( -2 );
125 }
126 if (UseDataInPlace_ != 1) {
127 assert ( RowMatrixA_ != 0 ) ;
128 if ( AddZeroToDiag_==0 && numentries_ != RowMatrixA_->NumGlobalNonzeros64()) {
129 std::cerr << " The number of non zero entries in the matrix has changed since the last call to SymbolicFactorization(). " ;
130 AMESOS_CHK_ERR( -3 );
131 }
132 assert ( ImportToSerial_.get() != 0 ) ;
134 *ImportToSerial_, Insert ));
135
136 if (AddZeroToDiag_ ) {
137 int OriginalTracebackMode = SerialCrsMatrixA_->GetTracebackMode() ;
138 SerialCrsMatrixA_->SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) ) ; // ExportToSerial is called both by PerformSymbolicFactorization() and PerformNumericFactorization(). When called by the latter, the call to insertglobalvalues is both unnecessary and illegal. Fortunately, Epetra allows us to ignore the error message by setting the traceback mode to 0.
139
140 //
141 // Add 0.0 to each diagonal entry to avoid empty diagonal entries;
142 //
143 double zero = 0.0;
144
145#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
146 if(SerialMap_->GlobalIndicesInt()) {
147 for ( int i = 0 ; i < SerialMap_->NumGlobalElements(); i++ ) {
148 if ( SerialCrsMatrixA_->LRID(i) >= 0 )
149 SerialCrsMatrixA_->InsertGlobalValues( i, 1, &zero, &i ) ;
150 }
151 }
152 else
153#endif
154#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
155 if(SerialMap_->GlobalIndicesLongLong()) {
156 for ( long long i = 0 ; i < SerialMap_->NumGlobalElements64(); i++ ) {
157 if ( SerialCrsMatrixA_->LRID(i) >= 0 )
158 SerialCrsMatrixA_->InsertGlobalValues( i, 1, &zero, &i ) ;
159 }
160 }
161 else
162#endif
163 throw "Amesos_Klu::ExportToSerial: ERROR, GlobalIndices type unknown.";
164
165 SerialCrsMatrixA_->SetTracebackMode( OriginalTracebackMode ) ;
166 }
167 AMESOS_CHK_ERR(SerialCrsMatrixA_->FillComplete());
168 AMESOS_CHK_ERR(SerialCrsMatrixA_->OptimizeStorage());
169
170 if( ! AddZeroToDiag_ && numentries_ != SerialMatrix_->NumGlobalNonzeros64()) {
171 std::cerr << " Amesos_Klu cannot handle this matrix. " ;
172 if ( Reindex_ ) {
173 std::cerr << "Unknown error" << std::endl ;
174 AMESOS_CHK_ERR( -4 );
175 } else {
176 std::cerr << " Try setting the Reindex parameter to true. " << std::endl ;
177#ifndef HAVE_AMESOS_EPETRAEXT
178 std::cerr << " You will need to rebuild the Amesos library with the EpetraExt library to use the reindex feature" << std::endl ;
179 std::cerr << " To rebuild Amesos with EpetraExt, add --enable-epetraext to your configure invocation" << std::endl ;
180#endif
181 AMESOS_CHK_ERR( -5 );
182 }
183 }
184 }
185
186 return 0;
187}
188//=============================================================================
189//
190// CreateLocalMatrixAndExporters() is called only by SymbolicFactorization()
191// for CrsMatrix objects. All objects should be created here. No assumptions
192// are made about the input operator. I.e. it can be completely different from
193// the matrix at the time of the previous call to SymbolicFactorization().
194//
195
197{
198 ResetTimer(0);
199
200 RowMatrixA_ = Problem_->GetMatrix(); // MS, 18-Apr-06 //
201 if (RowMatrixA_ == 0) AMESOS_CHK_ERR(-1);
202
203 const Epetra_Map &OriginalMatrixMap = RowMatrixA_->RowMatrixRowMap() ;
204 const Epetra_Map &OriginalDomainMap =
205 UseTranspose()?GetProblem()->GetOperator()->OperatorRangeMap():
206 GetProblem()->GetOperator()->OperatorDomainMap();
207 const Epetra_Map &OriginalRangeMap =
208 UseTranspose()?GetProblem()->GetOperator()->OperatorDomainMap():
209 GetProblem()->GetOperator()->OperatorRangeMap();
210
211 NumGlobalElements_ = RowMatrixA_->NumGlobalRows64();
212 numentries_ = RowMatrixA_->NumGlobalNonzeros64();
213 assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols64() );
214
215 //
216 // Create a serial matrix
217 //
218 assert( NumGlobalElements_ == OriginalMatrixMap.NumGlobalElements64() ) ;
219 int NumMyElements_ = 0 ;
220 if (MyPID_==0) NumMyElements_ = static_cast<int>(NumGlobalElements_);
221 //
222 // UseDataInPlace_ is set to 1 (true) only if everything is perfectly
223 // normal. Anything out of the ordinary reverts to the more expensive
224 // path.
225 //
226 UseDataInPlace_ = ( OriginalMatrixMap.NumMyElements() ==
227 OriginalMatrixMap.NumGlobalElements64() )?1:0;
228 if ( ! OriginalRangeMap.SameAs( OriginalMatrixMap ) ) UseDataInPlace_ = 0 ;
229 if ( ! OriginalDomainMap.SameAs( OriginalMatrixMap ) ) UseDataInPlace_ = 0 ;
230 if ( AddZeroToDiag_ ) UseDataInPlace_ = 0 ;
231 Comm().Broadcast( &UseDataInPlace_, 1, 0 ) ;
232
233 //
234 // Reindex matrix if necessary (and possible - i.e. CrsMatrix)
235 //
236 // For now, since I don't know how to determine if we need to reindex the matrix,
237 // we allow the user to control re-indexing through the "Reindex" parameter.
238 //
239 CrsMatrixA_ = dynamic_cast<Epetra_CrsMatrix *>(Problem_->GetOperator());
240
241 if ( Reindex_ ) {
242 if ( CrsMatrixA_ == 0 ) AMESOS_CHK_ERR(-7);
243#ifdef HAVE_AMESOS_EPETRAEXT
244 const Epetra_Map& OriginalMap = CrsMatrixA_->RowMap();
245 StdIndex_ = rcp( new Amesos_StandardIndex( OriginalMap ) );
246 //const Epetra_Map& OriginalColMap = CrsMatrixA_->RowMap();
247 StdIndexDomain_ = rcp( new Amesos_StandardIndex( OriginalDomainMap ) );
248 StdIndexRange_ = rcp( new Amesos_StandardIndex( OriginalRangeMap ) );
249
250 StdIndexMatrix_ = StdIndex_->StandardizeIndex( CrsMatrixA_ );
251#else
252 std::cerr << "Amesos_Klu requires EpetraExt to reindex matrices." << std::endl
253 << " Please rebuild with the EpetraExt library by adding --enable-epetraext to your configure invocation" << std::endl ;
254 AMESOS_CHK_ERR(-8);
255#endif
256 } else {
258 }
259
260 //
261 // At this point, StdIndexMatrix_ points to a matrix with
262 // standard indexing.
263 //
264
265 //
266 // Convert Original Matrix to Serial (if it is not already)
267 //
268 if (UseDataInPlace_ == 1) {
270 } else {
271#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
272 if(OriginalMatrixMap.GlobalIndicesInt()) {
273 int indexBase = OriginalMatrixMap.IndexBase();
274 SerialMap_ = rcp(new Epetra_Map((int) NumGlobalElements_, NumMyElements_, indexBase, Comm()));
275 }
276 else
277#endif
278#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
279 if(OriginalMatrixMap.GlobalIndicesLongLong()) {
280 long long indexBase = OriginalMatrixMap.IndexBase64();
281 SerialMap_ = rcp(new Epetra_Map((long long) NumGlobalElements_, NumMyElements_, indexBase, Comm()));
282 }
283 else
284#endif
285 throw "Amesos_Klu::CreateLocalMatrixAndExporters: Unknown Global Indices";
286
287 if ( Problem_->GetRHS() )
288 NumVectors_ = Problem_->GetRHS()->NumVectors() ;
289 else
290 NumVectors_ = 1 ;
291 SerialXextract_ = rcp( new Epetra_MultiVector(*SerialMap_,NumVectors_));
292 SerialBextract_ = rcp (new Epetra_MultiVector(*SerialMap_,NumVectors_));
293
294 ImportToSerial_ = rcp(new Epetra_Import ( *SerialMap_, StdIndexMatrix_->RowMatrixRowMap() ) );
295
296 if (ImportToSerial_.get() == 0) AMESOS_CHK_ERR(-9);
297
298 // Build the vector data import/export once and only once
299#define CHRIS
300#ifdef CHRIS
301 if(StdIndexMatrix_->RowMatrixRowMap().SameAs(StdIndexMatrix_->OperatorRangeMap()))
303 else
304 ImportRangeToSerial_ = rcp(new Epetra_Import ( *SerialMap_, StdIndexMatrix_->OperatorRangeMap() ) );
305
306 if(StdIndexMatrix_->RowMatrixRowMap().SameAs(StdIndexMatrix_->OperatorDomainMap()))
308 else
309 ImportDomainToSerial_ = rcp(new Epetra_Import ( *SerialMap_, StdIndexMatrix_->OperatorDomainMap() ) );
310#endif
311
312 SerialCrsMatrixA_ = rcp( new Epetra_CrsMatrix(Copy, *SerialMap_, 0) ) ;
314 }
315
316 MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_, 0);
317
318 return(0);
319}
320
321//=============================================================================
322//
323// See also pre and post conditions in Amesos_Klu.h
324// Preconditions:
325// firsttime specifies that this is the first time that
326// ConertToKluCrs has been called, i.e. in symbolic factorization.
327// No data allocation should happen unless firsttime=true.
328// SerialMatrix_ points to the matrix to be factored and solved
329// NumGlobalElements_ has been set to the dimension of the matrix
330// numentries_ has been set to the number of non-zeros in the matrix
331// (i.e. CreateLocalMatrixAndExporters() has been callded)
332//
333// Postconditions:
334// Ap, VecAi, VecAval contain the matrix as Klu needs it
335//
336//
338{
339 ResetTimer(0);
340
341 //
342 // Convert matrix to the form that Klu expects (Ap, VecAi, VecAval)
343 //
344
345 if (MyPID_==0) {
346 assert( NumGlobalElements_ == SerialMatrix_->NumGlobalRows64());
347 assert( NumGlobalElements_ == SerialMatrix_->NumGlobalCols64());
348 if ( ! AddZeroToDiag_ ) {
349 assert( numentries_ == SerialMatrix_->NumGlobalNonzeros64()) ;
350 } else {
351 numentries_ = SerialMatrix_->NumGlobalNonzeros64() ;
352 }
353
354 Epetra_CrsMatrix *CrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(SerialMatrix_);
355 bool StorageOptimized = ( CrsMatrix != 0 && CrsMatrix->StorageOptimized() );
356
357 if ( AddToDiag_ != 0.0 ) StorageOptimized = false ;
358
359 if ( firsttime ) {
360 Ap.resize( NumGlobalElements_+1 );
361 if ( ! StorageOptimized ) {
362 VecAi.resize( EPETRA_MAX( NumGlobalElements_, numentries_) ) ;
363 VecAval.resize( EPETRA_MAX( NumGlobalElements_, numentries_) ) ;
364 Ai = &VecAi[0];
365 Aval = &VecAval[0];
366 }
367 }
368
369 double *RowValues;
370 int *ColIndices;
371 int NumEntriesThisRow;
372
373 if( StorageOptimized ) {
374 if ( firsttime ) {
375 Ap[0] = 0;
376 for ( int MyRow = 0; MyRow <NumGlobalElements_; MyRow++ ) {
377 if( CrsMatrix->
378 ExtractMyRowView( MyRow, NumEntriesThisRow, RowValues,
379 ColIndices ) != 0 )
380 AMESOS_CHK_ERR( -10 );
381 if ( MyRow == 0 ) {
382 Ai = ColIndices ;
383 Aval = RowValues ;
384 }
385 Ap[MyRow+1] = Ap[MyRow] + NumEntriesThisRow ;
386 }
387 }
388 } else {
389
390 int Ai_index = 0 ;
391 int MyRow;
392
393 int MaxNumEntries_ = SerialMatrix_->MaxNumEntries();
394 if ( firsttime && CrsMatrix == 0 ) {
395 ColIndicesV_.resize(MaxNumEntries_);
396 RowValuesV_.resize(MaxNumEntries_);
397 }
398
399 for ( MyRow = 0; MyRow <NumGlobalElements_; MyRow++ ) {
400 if ( CrsMatrix != 0 ) {
401 if( CrsMatrix->
402 ExtractMyRowView( MyRow, NumEntriesThisRow, RowValues,
403 ColIndices ) != 0 )
404 AMESOS_CHK_ERR( -11 );
405
406 } else {
407 if( SerialMatrix_->
408 ExtractMyRowCopy( MyRow, MaxNumEntries_,
409 NumEntriesThisRow, &RowValuesV_[0],
410 &ColIndicesV_[0] ) != 0 )
411 AMESOS_CHK_ERR( -12 );
412
413 RowValues = &RowValuesV_[0];
414 ColIndices = &ColIndicesV_[0];
415 }
416
417 if ( firsttime ) {
418 Ap[MyRow] = Ai_index ;
419 for ( int j = 0; j < NumEntriesThisRow; j++ ) {
420 VecAi[Ai_index] = ColIndices[j] ;
421 // assert( VecAi[Ai_index] == Ai[Ai_index] ) ;
422 VecAval[Ai_index] = RowValues[j] ; // We have to do this because of the hacks to get aroun bug #1502
423 if (ColIndices[j] == MyRow) {
424 VecAval[Ai_index] += AddToDiag_;
425 }
426 Ai_index++;
427 }
428 } else {
429 for ( int j = 0; j < NumEntriesThisRow; j++ ) {
430 VecAval[Ai_index] = RowValues[j] ;
431 if (ColIndices[j] == MyRow) {
432 VecAval[Ai_index] += AddToDiag_;
433 }
434 Ai_index++;
435 }
436 }
437 }
438 Ap[MyRow] = Ai_index ;
439 }
440
441 }
442
443 MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
444
445 return 0;
446}
447
448//=============================================================================
449int Amesos_Klu::SetParameters( Teuchos::ParameterList &ParameterList ) {
450
451 // ========================================= //
452 // retrive KLU's parameters from list. //
453 // default values defined in the constructor //
454 // ========================================= //
455
456 // retrive general parameters
457
458 SetStatusParameters( ParameterList );
459 SetControlParameters( ParameterList );
460
461 if (ParameterList.isParameter("TrustMe") )
462 TrustMe_ = ParameterList.get<bool>( "TrustMe" );
463
464#if 0
465
466 unused for now
467
468 if (ParameterList.isSublist("Klu") ) {
469 Teuchos::ParameterList KluParams = ParameterList.sublist("Klu") ;
470 }
471#endif
472
473 return 0;
474}
475
476
477//=============================================================================
479{
480 ResetTimer(0);
481
482 int symbolic_ok = 0;
483
484 if (MyPID_ == 0) {
485
486 PrivateKluData_->common_ = rcp(new trilinos_klu_common());
487 trilinos_klu_defaults(&*PrivateKluData_->common_) ;
488 PrivateKluData_->Symbolic_ =
489 rcpWithDealloc(
490 trilinos_klu_analyze (static_cast<int>(NumGlobalElements_), &Ap[0], Ai, &*PrivateKluData_->common_ )
491 ,deallocFunctorDeleteWithCommon<trilinos_klu_symbolic>(PrivateKluData_->common_,trilinos_klu_free_symbolic)
492 ,true
493 );
494
495 symbolic_ok = (PrivateKluData_->Symbolic_.get() != NULL)
496 && PrivateKluData_->common_->status == TRILINOS_KLU_OK ;
497
498 }
499
500 // Communicate the state of the symbolic factorization with everyone.
501 Comm().Broadcast(&symbolic_ok, 1, 0);
502
503 if ( !symbolic_ok ) {
504 if (MyPID_ == 0) {
506 } else
508 }
509
510 SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_, 0);
511
512 return 0;
513}
514
515//=============================================================================
517{
518
519 // Changed this; it was "if (!TrustMe)...
520 // The behavior is not intuitive. Maybe we should introduce a new
521 // parameter, FastSolvers or something like that, that does not perform
522 // any AddTime, ResetTimer, GetTime.
523 // HKT, 1/18/2007: The "TrustMe_" flag was put back in this code due to performance degradation; Bug# 3042.
524
525 if (! TrustMe_ ) ResetTimer(0);
526
527 // This needs to be 1 just in case the pivot ordering is reused.
528 int numeric_ok = 1;
529
530 if (MyPID_ == 0) {
531
532 bool factor_with_pivoting = true ;
533
534 // set the default parameters
535 PrivateKluData_->common_->scale = ScaleMethod_ ;
536
537 const bool NumericNonZero = PrivateKluData_->Numeric_.get() != 0 ;
538
539 // see if we can "refactorize"
540 if ( refactorize_ && NumericNonZero ) {
541 // refactorize using the existing Symbolic and Numeric objects, and
542 // using the identical pivot ordering as the prior klu_factor.
543 // No partial pivoting is done.
544 int result = trilinos_klu_refactor (&Ap[0], Ai, Aval,
545 &*PrivateKluData_->Symbolic_,
546 &*PrivateKluData_->Numeric_, &*PrivateKluData_->common_) ;
547 // Did it work?
548 const bool refactor_ok = result == 1 && PrivateKluData_->common_->status == TRILINOS_KLU_OK ;
549 if ( refactor_ok ) {
550
551 trilinos_klu_rcond (&*PrivateKluData_->Symbolic_,
552 &*PrivateKluData_->Numeric_,
553 &*PrivateKluData_->common_) ;
554
555 double rcond = PrivateKluData_->common_->rcond;
556
557 if ( rcond > rcond_threshold_ ) {
558 // factorizing without pivot worked fine. We are done.
559 factor_with_pivoting = false ;
560 }
561 }
562 }
563
564 if ( factor_with_pivoting ) {
565
566 // factor with partial pivoting:
567 // either this is the first time we are factoring the matrix, or the
568 // refactorize parameter is false, or we tried to refactorize and
569 // found it to be too inaccurate.
570
571 // factor the matrix using partial pivoting
572 PrivateKluData_->Numeric_ =
573 rcpWithDealloc( trilinos_klu_factor(&Ap[0], Ai, Aval,
574 &*PrivateKluData_->Symbolic_, &*PrivateKluData_->common_),
575 deallocFunctorDeleteWithCommon<trilinos_klu_numeric>(PrivateKluData_->common_,trilinos_klu_free_numeric)
576 ,true
577 );
578
579 numeric_ok = PrivateKluData_->Numeric_.get()!=NULL
580 && PrivateKluData_->common_->status == TRILINOS_KLU_OK ;
581
582 }
583 }
584
585 // Communicate the state of the numeric factorization with everyone.
586 Comm().Broadcast(&numeric_ok, 1, 0);
587
588 if ( ! numeric_ok ) {
589 if (MyPID_ == 0) {
591 } else
593 }
594
595 if ( !TrustMe_ ) {
596 NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
597 }
598
599 return 0;
600}
601
602//=============================================================================
604 bool OK = true;
605
606 // Comment by Tim: The following code seems suspect. The variable "OK"
607 // is not set if the condition is true.
608 // Does the variable "OK" default to true?
609 if ( GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints64() !=
610 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints64() ) {
611 OK = false;
612 }
613 return OK;
614}
615
616//=============================================================================
618{
619 MyPID_ = Comm().MyPID();
620 NumProcs_ = Comm().NumProc();
621
624
625 CreateTimer(Comm(), 2);
626
627 ResetTimer(1);
628
629 RowMatrixA_ = Problem_->GetMatrix();
630 if (RowMatrixA_->NumGlobalRows64() == 0 || RowMatrixA_->NumGlobalCols64()
631 == 0)
632 {
635 return 0;
636 }
637
638 // "overhead" time for the following method is considered here
640 assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols64() );
641
642 //
643 // Perform checks in SymbolicFactorization(), but none in
644 // NumericFactorization() or Solve()
645 //
646 if ( TrustMe_ ) {
647 if ( CrsMatrixA_ == 0 ) AMESOS_CHK_ERR( -15 );
648 if( UseDataInPlace_ != 1 ) AMESOS_CHK_ERR( -15 ) ;
649 if( Reindex_ ) AMESOS_CHK_ERR( -15 ) ;
650 if( ! Problem_->GetLHS() ) AMESOS_CHK_ERR( -15 ) ;
651 if( ! Problem_->GetRHS() ) AMESOS_CHK_ERR( -15 ) ;
652 if( ! Problem_->GetLHS()->NumVectors() ) AMESOS_CHK_ERR( -15 ) ;
653 if( ! Problem_->GetRHS()->NumVectors() ) AMESOS_CHK_ERR( -15 ) ;
654 SerialB_ = Teuchos::rcp(Problem_->GetRHS(),false) ;
655 SerialX_ = Teuchos::rcp(Problem_->GetLHS(),false) ;
656 NumVectors_ = SerialX_->NumVectors();
657 if (MyPID_ == 0) {
660 }
661 }
662
663 // "overhead" time for the following two methods is considered here
665
667
668 OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
669
670 // All this time if KLU time
672
674
676
677 return 0;
678}
679
680//=============================================================================
682{
683 RowMatrixA_ = Problem_->GetMatrix();
684 if (RowMatrixA_->NumGlobalRows64() == 0 || RowMatrixA_->NumGlobalCols64()
685 == 0)
686 {
689 return 0;
690 }
691
692 if ( !TrustMe_ )
693 {
695 if (IsSymbolicFactorizationOK_ == false)
697
698 ResetTimer(1); // "overhead" time
699
700 Epetra_CrsMatrix *CrsMatrixA = dynamic_cast<Epetra_CrsMatrix *>(RowMatrixA_);
701 if ( CrsMatrixA == 0 ) // hack to get around Bug #1502
703 assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols64() );
704 if ( AddZeroToDiag_ == 0 )
705 assert( numentries_ == RowMatrixA_->NumGlobalNonzeros64() );
706
708
709 if ( CrsMatrixA == 0 ) { // continuation of hack to avoid bug #1502
711 } else {
713 }
714
715 OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
716 }
717
718 // this time is all for KLU
720
722
724
725 return 0;
726}
727
728//=============================================================================
730{
731 RowMatrixA_ = Problem_->GetMatrix();
732 if (RowMatrixA_->NumGlobalRows64() == 0 || RowMatrixA_->NumGlobalCols64()
733 == 0)
734 {
735 ++NumSolve_;
736 return 0;
737 }
738
739 Epetra_MultiVector* vecX = 0 ;
740 Epetra_MultiVector* vecB = 0 ;
741
742#ifdef HAVE_AMESOS_EPETRAEXT
743 Teuchos::RCP<Epetra_MultiVector> vecX_rcp;
744 Teuchos::RCP<Epetra_MultiVector> vecB_rcp;
745#endif
746
747#ifdef Bug_8212
748 // This demonstrates Bug #2812 - Valgrind does not catch this
749 // memory leak
750 lose_this_ = (int *) malloc( 300 ) ;
751
752#ifdef Bug_8212_B
753 // This demonstrates Bug #2812 - Valgrind does catch this
754 // use of unitialized data - but only in TestOptions/TestOptions.exe
755 // not in Test_Basic/amesos_test.exe
756 //
757 if ( lose_this_[0] == 12834 ) {
758 std::cout << __FILE__ << "::" << __LINE__ << " very unlikely to happen " << std::endl ;
759 }
760#endif
761#endif
762
763 if ( !TrustMe_ ) {
764
765 SerialB_ = Teuchos::rcp(Problem_->GetRHS(),false);
766 SerialX_ = Teuchos::rcp(Problem_->GetLHS(),false);
767
768 Epetra_MultiVector* OrigVecX ;
769 Epetra_MultiVector* OrigVecB ;
770
771 if (IsNumericFactorizationOK_ == false)
773
774 ResetTimer(1);
775
776 //
777 // Reindex the LHS and RHS
778 //
779 OrigVecX = Problem_->GetLHS() ;
780 OrigVecB = Problem_->GetRHS() ;
781
782 if ( Reindex_ ) {
783#ifdef HAVE_AMESOS_EPETRAEXT
784 vecX_rcp = StdIndexDomain_->StandardizeIndex( *OrigVecX ) ;
785 vecB_rcp = StdIndexRange_->StandardizeIndex( *OrigVecB ) ;
786
787 vecX = &*vecX_rcp;
788 vecB = &*vecB_rcp;
789#else
790 AMESOS_CHK_ERR( -13 ) ; // Amesos_Klu can't handle non-standard indexing without EpetraExt
791#endif
792 } else {
793 vecX = OrigVecX ;
794 vecB = OrigVecB ;
795 }
796
797 if ((vecX == 0) || (vecB == 0))
798 AMESOS_CHK_ERR(-1); // something wrong in input
799
800 // Extract Serial versions of X and B
801
802 ResetTimer(0);
803
804 // Copy B to the serial version of B
805 //
806 if (UseDataInPlace_ == 1) {
807#ifdef HAVE_AMESOS_EPETRAEXT
808 if(vecX_rcp==Teuchos::null)
809 SerialX_ = Teuchos::rcp(vecX,false);
810 else
811 SerialX_ = vecX_rcp;
812
813 if(vecB_rcp==Teuchos::null)
814 SerialB_ = Teuchos::rcp(vecB,false);
815 else
816 SerialB_ = vecB_rcp;
817#else
818 SerialB_ = Teuchos::rcp(vecB,false);
819 SerialX_ = Teuchos::rcp(vecX,false);
820#endif
821 NumVectors_ = Problem_->GetRHS()->NumVectors() ;
822 } else {
823 assert (UseDataInPlace_ == 0);
824
825 if( NumVectors_ != Problem_->GetRHS()->NumVectors() ) {
826 NumVectors_ = Problem_->GetRHS()->NumVectors() ;
827 SerialXextract_ = rcp( new Epetra_MultiVector(*SerialMap_,NumVectors_));
828 SerialBextract_ = rcp (new Epetra_MultiVector(*SerialMap_,NumVectors_));
829 }
830 if (NumVectors_ != vecB->NumVectors())
831 AMESOS_CHK_ERR(-1); // internal error
832
833 //ImportRangeToSerial_ = rcp(new Epetra_Import ( *SerialMap_, vecB->Map() ) );
834 //if ( SerialBextract_->Import(*vecB,*ImportRangeToSerial_,Insert) )
835 Epetra_Import *UseImport;
836 if(!UseTranspose_) UseImport=&*ImportRangeToSerial_;
837 else UseImport=&*ImportDomainToSerial_;
838 if ( SerialBextract_->Import(*vecB,*UseImport,Insert) )
839 AMESOS_CHK_ERR( -1 ) ; // internal error
840
841 SerialB_ = Teuchos::rcp(&*SerialBextract_,false) ;
842 SerialX_ = Teuchos::rcp(&*SerialXextract_,false) ;
843 }
844
845 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
846
847 // Call KLU to perform the solve
848
849 ResetTimer(0);
850 if (MyPID_ == 0) {
854 AMESOS_CHK_ERR(-1);
855 }
856
857 OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
858 }
859 else
860 {
861 SerialB_ = Teuchos::rcp(Problem_->GetRHS(),false) ;
862 SerialX_ = Teuchos::rcp(Problem_->GetLHS(),false) ;
863 NumVectors_ = SerialX_->NumVectors();
864 if (MyPID_ == 0) {
867 }
868 }
869
870 if ( MyPID_ == 0) {
871 if ( NumVectors_ == 1 ) {
872 for ( int i = 0 ; i < NumGlobalElements_ ; i++ )
874 } else {
875 SerialX_->Scale(1.0, *SerialB_ ) ; // X = B (Klu overwrites B with X)
876 }
877 if (UseTranspose()) {
878 trilinos_klu_solve( &*PrivateKluData_->Symbolic_, &*PrivateKluData_->Numeric_,
880 } else {
881 trilinos_klu_tsolve( &*PrivateKluData_->Symbolic_, &*PrivateKluData_->Numeric_,
883 }
884 }
885
886 if ( !TrustMe_ ) {
887 SolveTime_ = AddTime("Total solve time", SolveTime_, 0);
888
889 // Copy X back to the original vector
890
891 ResetTimer(0);
892 ResetTimer(1);
893
894 if (UseDataInPlace_ == 0) {
895 Epetra_Import *UseImport;
896 if(!UseTranspose_) UseImport=&*ImportDomainToSerial_;
897 else UseImport=&*ImportRangeToSerial_;
898 // ImportDomainToSerial_ = rcp(new Epetra_Import ( *SerialMap_, vecX->Map() ) );
899 vecX->Export( *SerialX_, *UseImport, Insert ) ;
900
901 } // otherwise we are already in place.
902
903 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
904
905#if 0
906 //
907 // ComputeTrueResidual causes TestOptions to fail on my linux box
908 // Bug #1417
910 ComputeTrueResidual(*SerialMatrix_, *vecX, *vecB, UseTranspose(), "Amesos_Klu");
911#endif
912
914 ComputeVectorNorms(*vecX, *vecB, "Amesos_Klu");
915
916 OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
917
918 }
919 ++NumSolve_;
920
921 return(0) ;
922}
923
924// ================================================ ====== ==== ==== == =
925
927{
928
929 if (MyPID_) return;
930
931 PrintLine();
932
933 std::cout << "Amesos_Klu : Matrix has " << NumGlobalElements_ << " rows"
934 << " and " << numentries_ << " nonzeros" << std::endl;
935 std::cout << "Amesos_Klu : Nonzero elements per row = "
936 << 1.0*numentries_/NumGlobalElements_ << std::endl;
937 std::cout << "Amesos_Klu : Percentage of nonzero elements = "
938 << 100.0*numentries_/(pow(double(NumGlobalElements_),double(2.0))) << std::endl;
939 std::cout << "Amesos_Klu : Use transpose = " << UseTranspose_ << std::endl;
940
941 PrintLine();
942
943 return;
944
945}
946
947// ================================================ ====== ==== ==== == =
948
950{
951 if (MyPID_) return;
952
953 double ConTime = GetTime(MtxConvTime_);
954 double MatTime = GetTime(MtxRedistTime_);
955 double VecTime = GetTime(VecRedistTime_);
956 double SymTime = GetTime(SymFactTime_);
957 double NumTime = GetTime(NumFactTime_);
958 double SolTime = GetTime(SolveTime_);
959 double OveTime = GetTime(OverheadTime_);
960
962 SymTime /= NumSymbolicFact_;
963
964 if (NumNumericFact_)
965 NumTime /= NumNumericFact_;
966
967 if (NumSolve_)
968 SolTime /= NumSolve_;
969
970 std::string p = "Amesos_Klu : ";
971 PrintLine();
972
973 std::cout << p << "Time to convert matrix to Klu format = "
974 << ConTime << " (s)" << std::endl;
975 std::cout << p << "Time to redistribute matrix = "
976 << MatTime << " (s)" << std::endl;
977 std::cout << p << "Time to redistribute vectors = "
978 << VecTime << " (s)" << std::endl;
979 std::cout << p << "Number of symbolic factorizations = "
980 << NumSymbolicFact_ << std::endl;
981 std::cout << p << "Time for sym fact = "
982 << SymTime * NumSymbolicFact_ << " (s), avg = " << SymTime << " (s)" << std::endl;
983 std::cout << p << "Number of numeric factorizations = "
984 << NumNumericFact_ << std::endl;
985 std::cout << p << "Time for num fact = "
986 << NumTime * NumNumericFact_ << " (s), avg = " << NumTime << " (s)" << std::endl;
987 std::cout << p << "Number of solve phases = "
988 << NumSolve_ << std::endl;
989 std::cout << p << "Time for solve = "
990 << SolTime * NumSolve_ << " (s), avg = " << SolTime << " (s)" << std::endl;
991
992 double tt = SymTime * NumSymbolicFact_ + NumTime * NumNumericFact_ + SolTime * NumSolve_;
993 if (tt != 0)
994 {
995 std::cout << p << "Total time spent in Amesos = " << tt << " (s) " << std::endl;
996 std::cout << p << "Total time spent in the Amesos interface = " << OveTime << " (s)" << std::endl;
997 std::cout << p << "(the above time does not include KLU time)" << std::endl;
998 std::cout << p << "Amesos interface time / total time = " << OveTime / tt << std::endl;
999 }
1000
1001 PrintLine();
1002
1003 return;
1004}
const int NumericallySingularMatrixError
const int StructurallySingularMatrixError
#define AMESOS_CHK_ERR(a)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
double rcond_threshold_
If error is greater than this value, perform symbolic and numeric factorization with full partial piv...
bool Reindex_
If true, the Amesos class should reindex the matrix to standard indexing (i.e.
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.
Teuchos::RCP< trilinos_klu_symbolic > Symbolic_
Teuchos::RCP< trilinos_klu_numeric > Numeric_
Teuchos::RCP< trilinos_klu_common > common_
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
void PrintTiming() const
Prints timing information.
const Epetra_LinearProblem * GetProblem() const
Get a pointer to the Problem.
Definition Amesos_Klu.h:147
Teuchos::RCP< Epetra_MultiVector > SerialXextract_
Serial versions of the LHS and RHS (if necessary)
Definition Amesos_Klu.h:324
double * Aval
Definition Amesos_Klu.h:279
Teuchos::RCP< Epetra_MultiVector > SerialB_
Serial versions of the LHS and RHS (may point to the original vector if serial)
Definition Amesos_Klu.h:321
int MtxConvTime_
Definition Amesos_Klu.h:341
int PerformNumericFactorization()
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int SerialXlda_
Definition Amesos_Klu.h:258
std::vector< double > VecAval
Definition Amesos_Klu.h:278
int CreateLocalMatrixAndExporters()
bool TrustMe_
If true, no checks are made and the matrix is assume to be distributed.
Definition Amesos_Klu.h:314
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Definition Amesos_Klu.h:164
bool UseTranspose() const
Returns the current UseTranspose setting.
Definition Amesos_Klu.h:162
bool UseTranspose_
If true, the transpose of A is used.
Definition Amesos_Klu.h:328
Epetra_CrsMatrix * CrsMatrixA_
Operator converted to a CrsMatrix.
Definition Amesos_Klu.h:292
~Amesos_Klu(void)
Amesos_Klu Destructor.
Teuchos::RCP< Epetra_Import > ImportDomainToSerial_
Definition Amesos_Klu.h:339
Teuchos::RCP< Epetra_Import > ImportToSerial_
Importer to process 0.
Definition Amesos_Klu.h:337
Teuchos::RCP< Epetra_Map > SerialMap_
Points to a Serial Map (unused if UseDataInPlace_ == 1 )
Definition Amesos_Klu.h:302
std::vector< int > ColIndicesV_
Only used for RowMatrices to extract copies.
Definition Amesos_Klu.h:333
Teuchos::RCP< Amesos_Klu_Pimpl > PrivateKluData_
Definition Amesos_Klu.h:267
int MtxRedistTime_
Quick access ids for the individual timings.
Definition Amesos_Klu.h:341
int ExportToSerial()
Teuchos::RCP< Amesos_StandardIndex > StdIndexDomain_
Definition Amesos_Klu.h:270
Teuchos::RCP< Epetra_MultiVector > SerialX_
Definition Amesos_Klu.h:322
double * SerialBvalues_
Definition Amesos_Klu.h:319
int SymFactTime_
Definition Amesos_Klu.h:342
int UseDataInPlace_
1 if Problem_->GetOperator() is stored entirely on process 0
Definition Amesos_Klu.h:283
void PrintStatus() const
Prints information about the factorization and solution phases.
std::vector< double > RowValuesV_
Only used for RowMatrices to extract copies.
Definition Amesos_Klu.h:335
Epetra_RowMatrix * RowMatrixA_
Operator converted to a RowMatrix.
Definition Amesos_Klu.h:290
Amesos_Klu(const Epetra_LinearProblem &LinearProblem)
Amesos_Klu Constructor.
Teuchos::RCP< Epetra_Import > ImportRangeToSerial_
Definition Amesos_Klu.h:338
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrixA_
Points to a Serial Copy of A (unused if UseDataInPlace_==1)
Definition Amesos_Klu.h:304
Teuchos::RCP< Amesos_StandardIndex > StdIndex_
Definition Amesos_Klu.h:268
const Epetra_LinearProblem * Problem_
Pointer to the linear system problem.
Definition Amesos_Klu.h:330
Teuchos::RCP< Epetra_MultiVector > SerialBextract_
Definition Amesos_Klu.h:325
int NumFactTime_
Definition Amesos_Klu.h:342
std::vector< int > VecAi
Definition Amesos_Klu.h:277
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int VecRedistTime_
Definition Amesos_Klu.h:341
Teuchos::RCP< Amesos_StandardIndex > StdIndexRange_
Definition Amesos_Klu.h:269
int Solve()
Solves A X = B (or AT x = B)
int ConvertToKluCRS(bool firsttime)
int PerformSymbolicFactorization()
long long numentries_
Number of non-zero entries in Problem_->GetOperator()
Definition Amesos_Klu.h:285
std::vector< int > Ap
Ap, Ai, Aval form the compressed row storage used by Klu Ai and Aval can point directly into a matrix...
Definition Amesos_Klu.h:276
Epetra_RowMatrix * StdIndexMatrix_
Points to a Contiguous Copy of A.
Definition Amesos_Klu.h:306
bool MatrixShapeOK() const
Returns true if KLU can handle this matrix shape.
int NumVectors_
Number of vectors in RHS and LHS.
Definition Amesos_Klu.h:316
int OverheadTime_
Definition Amesos_Klu.h:342
Epetra_RowMatrix * SerialMatrix_
Points to a Serial Copy of A.
Definition Amesos_Klu.h:308
double * SerialXBvalues_
Pointer to the actual values in the serial version of X and B.
Definition Amesos_Klu.h:318
long long NumGlobalElements_
Number of rows and columns in the Problem_->GetOperator()
Definition Amesos_Klu.h:287
bool PrintTiming_
If true, prints timing information in the destructor.
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
int verbose_
Toggles the output level.
int NumSymbolicFact_
Number of symbolic factorization phases.
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
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)
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
int NumNumericFact_
Number of numeric factorization phases.
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 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.