Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Amesos_Paraklete.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_Paraklete.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"
36extern "C" {
37#include "amesos_paraklete_decl.h"
38}
39
40// The following hack allows assert to work even though it has been turned off by paraklete.h bug #1952
41
42#undef _ASSERT_H
43#undef assert
44#undef __ASSERT_VOID_CAST
45#undef NDEBUG
46#include <assert.h>
47
48namespace {
49
50using Teuchos::RCP;
51
52template<class T, class DeleteFunctor>
53class DeallocFunctorDeleteWithCommon
54{
55public:
56 DeallocFunctorDeleteWithCommon(
57 const RCP<paraklete_common> &common
58 ,DeleteFunctor deleteFunctor
59 )
60 : common_(common), deleteFunctor_(deleteFunctor)
61 {}
62 typedef T ptr_t;
63 void free( T* ptr ) {
64 if(ptr) deleteFunctor_(&ptr,&*common_);
65 }
66private:
67 Teuchos::RCP<paraklete_common> common_;
68 DeleteFunctor deleteFunctor_;
69 DeallocFunctorDeleteWithCommon(); // Not defined and not to be called!
70};
71
72template<class T, class DeleteFunctor>
73DeallocFunctorDeleteWithCommon<T,DeleteFunctor>
74deallocFunctorDeleteWithCommon(
75 const RCP<paraklete_common> &common
76 ,DeleteFunctor deleteFunctor
77 )
78{
79 return DeallocFunctorDeleteWithCommon<T,DeleteFunctor>(common,deleteFunctor);
80}
81
82
83} // namespace
84
85
86
88public:
89
90
91 cholmod_sparse pk_A_ ;
92 Teuchos::RCP<paraklete_symbolic> LUsymbolic_ ;
93 Teuchos::RCP<paraklete_numeric> LUnumeric_ ;
94 Teuchos::RCP<paraklete_common> common_;
95
96
97
98} ;
99
100//=============================================================================
101Amesos_Paraklete::Amesos_Paraklete(const Epetra_LinearProblem &prob ) :
102 PrivateParakleteData_( rcp( new Amesos_Paraklete_Pimpl() ) ),
103 ParakleteComm_(0),
104 CrsMatrixA_(0),
105 TrustMe_(false),
106 UseTranspose_(false),
107 Problem_(&prob),
108 MtxConvTime_(-1),
109 MtxRedistTime_(-1),
110 VecRedistTime_(-1),
111 SymFactTime_(-1),
112 NumFactTime_(-1),
113 SolveTime_(-1),
114 OverheadTime_(-1)
115{
116
117 // MS // move declaration of Problem_ above because I need it
118 // MS // set up before calling Comm()
119 MaxProcesses_ = -3; // Use all processes unless the user requests otherwise
120 Teuchos::ParameterList ParamList ;
121 SetParameters( ParamList ) ;
122}
123
124//=============================================================================
126
127
128 if(ParakleteComm_) {
129 MPI_Comm_free(&ParakleteComm_);
130 ParakleteComm_ = 0 ;
131 }
132
133 // print out some information if required by the user
134 if( (verbose_ && PrintTiming_) || verbose_ == 2 ) PrintTiming();
135 if( (verbose_ && PrintStatus_) || verbose_ == 2 ) PrintStatus();
136}
137
138//=============================================================================
140{
141 if ( AddZeroToDiag_==0 && numentries_ != RowMatrixA_->NumGlobalNonzeros()) {
142 std::cerr << " The number of non zero entries in the matrix has changed since the last call to SymbolicFactorization(). " ;
143 AMESOS_CHK_ERR( -2 );
144 }
145 if ( UseDataInPlace_ == 0 ) {
146 assert ( RowMatrixA_ != 0 ) ;
147 if ( AddZeroToDiag_==0 && numentries_ != RowMatrixA_->NumGlobalNonzeros()) {
148 std::cerr << " The number of non zero entries in the matrix has changed since the last call to SymbolicFactorization(). " ;
149 AMESOS_CHK_ERR( -2 );
150 }
151
152#ifdef HAVE_AMESOS_EPETRAEXT
153 if ( transposer_.get() != 0 ) {
154 // int OriginalTracebackMode = CrsMatrixA_->GetTracebackMode();
155 // CrsMatrixA_->SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) );
156 transposer_->fwd() ;
157 // CrsMatrixA_->SetTracebackMode( OriginalTracebackMode );
158 }
159#endif
160 assert ( ImportToSerial_.get() != 0 ) ;
162 *ImportToSerial_, Insert ));
163
164 if (AddZeroToDiag_ ) {
165 int OriginalTracebackMode = SerialCrsMatrixA_->GetTracebackMode() ;
166 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.
167
168 //
169 // Add 0.0 to each diagonal entry to avoid empty diagonal entries;
170 //
171 double zero = 0.0;
172 for ( int i = 0 ; i < SerialMap_->NumGlobalElements(); i++ )
173 if ( SerialCrsMatrixA_->LRID(i) >= 0 )
174 SerialCrsMatrixA_->InsertGlobalValues( i, 1, &zero, &i ) ;
175 SerialCrsMatrixA_->SetTracebackMode( OriginalTracebackMode ) ;
176 }
177 AMESOS_CHK_ERR(SerialCrsMatrixA_->FillComplete());
178 //AMESOS_CHK_ERR(SerialCrsMatrixA_->OptimizeStorage());
179
180 if( !AddZeroToDiag_ && numentries_ != SerialMatrix_->NumGlobalNonzeros()) {
181 std::cerr << " Amesos_Paraklete cannot handle this matrix. " ;
182 if ( Reindex_ ) {
183 std::cerr << "Unknown error" << std::endl ;
184 AMESOS_CHK_ERR( -5 );
185 } else {
186 std::cerr << " Try setting the Reindex parameter to true. " << std::endl ;
187#ifndef HAVE_AMESOS_EPETRAEXT
188 std::cerr << " You will need to rebuild the Amesos library with the EpetraExt library to use the reindex feature" << std::endl ;
189 std::cerr << " To rebuild Amesos with EpetraExt, add --enable-epetraext to your configure invocation" << std::endl ;
190#endif
191 AMESOS_CHK_ERR( -3 );
192 }
193 }
194 }
195 return 0;
196}
197//=============================================================================
198//
199// CreateLocalMatrixAndExporters() is called only by SymbolicFactorization()
200// for CrsMatrix objects. All objects should be created here. No assumptions
201// are made about the input operator. I.e. it can be completely different from
202// the matrix at the time of the previous call to SymbolicFactorization().
203//
204
206{
207 ResetTimer(0);
208
209 RowMatrixA_ = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
210 if (RowMatrixA_ == 0) AMESOS_CHK_ERR(-1);
211
212 const Epetra_Map &OriginalMatrixMap = RowMatrixA_->RowMatrixRowMap() ;
213 const Epetra_Map &OriginalDomainMap =
214 UseTranspose()?GetProblem()->GetOperator()->OperatorRangeMap():
215 GetProblem()->GetOperator()->OperatorDomainMap();
216 const Epetra_Map &OriginalRangeMap =
217 UseTranspose()?GetProblem()->GetOperator()->OperatorDomainMap():
218 GetProblem()->GetOperator()->OperatorRangeMap();
219
220 NumGlobalElements_ = RowMatrixA_->NumGlobalRows();
221 numentries_ = RowMatrixA_->NumGlobalNonzeros();
222 assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols() );
223
224 //
225 // Create a serial matrix
226 //
227 assert( NumGlobalElements_ == OriginalMatrixMap.NumGlobalElements() ) ;
228 int NumMyElements_ = 0 ;
229 if (MyPID_==0) NumMyElements_ = NumGlobalElements_;
230
231 //
232 // UseDataInPlace_ is set to 1 (true) only if everything is perfectly
233 // normal. Anything out of the ordinary reverts to the more expensive
234 // path.
235 //
236 UseDataInPlace_ = ( OriginalMatrixMap.NumMyElements() ==
237 OriginalMatrixMap.NumGlobalElements() )?1:0;
238 if ( ! OriginalRangeMap.SameAs( OriginalMatrixMap ) ) UseDataInPlace_ = 0 ;
239 if ( ! OriginalDomainMap.SameAs( OriginalMatrixMap ) ) UseDataInPlace_ = 0 ;
240 if ( AddZeroToDiag_ ) UseDataInPlace_ = 0 ;
241 Comm().Broadcast( &UseDataInPlace_, 1, 0 ) ;
242
243 UseDataInPlace_ = 0 ; // bug - remove this someday.
244
245 //
246 // Reindex matrix if necessary (and possible - i.e. CrsMatrix)
247 //
248 // For now, since I don't know how to determine if we need to reindex the matrix,
249 // we allow the user to control re-indexing through the "Reindex" parameter.
250 //
251 CrsMatrixA_ = dynamic_cast<Epetra_CrsMatrix *>(Problem_->GetOperator());
252 Epetra_CrsMatrix* CcsMatrixA = 0 ;
253
254 //
255 // CcsMatrixA points to a matrix in Compressed Column Format
256 // i.e. the format needed by Paraklete. If we are solving
257 // A^T x = b, CcsMatrixA = CrsMatrixA_, otherwise we must
258 // transpose the matrix.
259 //
260 if (UseTranspose()) {
261 CcsMatrixA = CrsMatrixA_ ;
262 } else {
263 if ( CrsMatrixA_ == 0 ) AMESOS_CHK_ERR( -7 ); // Amesos_Paraklete only supports CrsMatrix objects in the non-transpose case
264#ifdef HAVE_AMESOS_EPETRAEXT
265 bool MakeDataContiguous = true;
266 transposer_ = rcp ( new EpetraExt::RowMatrix_Transpose( MakeDataContiguous ));
267
268 int OriginalTracebackMode = CrsMatrixA_->GetTracebackMode();
269 CrsMatrixA_->SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) );
270
271 CcsMatrixA = &(dynamic_cast<Epetra_CrsMatrix&>(((*transposer_)(*CrsMatrixA_))));
272 CrsMatrixA_->SetTracebackMode( OriginalTracebackMode );
273#else
274 std::cerr << "Amesos_Paraklete requires the EpetraExt library to solve non-transposed problems. " << std::endl ;
275 std::cerr << " To rebuild Amesos with EpetraExt, add --enable-epetraext to your configure invocation" << std::endl ;
276 AMESOS_CHK_ERR( -3 );
277#endif
278 }
279
280
281 if ( Reindex_ ) {
282 if ( CcsMatrixA == 0 ) AMESOS_CHK_ERR(-4);
283#ifdef HAVE_AMESOS_EPETRAEXT
284 const Epetra_Map& OriginalMap = CcsMatrixA->RowMap();
285 StdIndex_ = rcp( new Amesos_StandardIndex( OriginalMap ) );
286 //const Epetra_Map& OriginalColMap = CcsMatrixA->RowMap();
287 StdIndexDomain_ = rcp( new Amesos_StandardIndex( OriginalDomainMap ) );
288 StdIndexRange_ = rcp( new Amesos_StandardIndex( OriginalRangeMap ) );
289
290 StdIndexMatrix_ = StdIndex_->StandardizeIndex( CcsMatrixA );
291#else
292 std::cerr << "Amesos_Paraklete requires EpetraExt to reindex matrices." << std::endl
293 << " Please rebuild with the EpetraExt library by adding --enable-epetraext to your configure invocation" << std::endl ;
294 AMESOS_CHK_ERR(-4);
295#endif
296 } else {
297 if ( UseTranspose() ){
299 } else {
300 StdIndexMatrix_ = CcsMatrixA ;
301 }
302 }
303
304 //
305 // At this point, StdIndexMatrix_ points to a matrix with
306 // standard indexing.
307 //
308
309 //
310 // Convert Original Matrix to Serial (if it is not already)
311 //
312 if (UseDataInPlace_ == 1) {
314 } else {
315 SerialMap_ = rcp(new Epetra_Map(NumGlobalElements_, NumMyElements_, 0, Comm()));
316
317 if ( Problem_->GetRHS() )
318 NumVectors_ = Problem_->GetRHS()->NumVectors() ;
319 else
320 NumVectors_ = 1 ;
321 SerialXextract_ = rcp( new Epetra_MultiVector(*SerialMap_,NumVectors_));
322 SerialBextract_ = rcp (new Epetra_MultiVector(*SerialMap_,NumVectors_));
323
324 ImportToSerial_ = rcp(new Epetra_Import ( *SerialMap_, StdIndexMatrix_->RowMatrixRowMap() ) );
325 if (ImportToSerial_.get() == 0) AMESOS_CHK_ERR(-5);
326
327 SerialCrsMatrixA_ = rcp( new Epetra_CrsMatrix(Copy, *SerialMap_, 0) ) ;
329 }
330
331 MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_, 0);
332
333 return(0);
334}
335
336//=============================================================================
337//
338// See also pre and post conditions in Amesos_Paraklete.h
339// Preconditions:
340// firsttime specifies that this is the first time that
341// ConertToParakleteCrs has been called, i.e. in symbolic factorization.
342// No data allocation should happen unless firsttime=true.
343// SerialMatrix_ points to the matrix to be factored and solved
344// NumGlobalElements_ has been set to the dimension of the matrix
345// numentries_ has been set to the number of non-zeros in the matrix
346// (i.e. CreateLocalMatrixAndExporters() has been callded)
347//
348// Postconditions:
349// pk_A_ contains the matrix as Paraklete needs it
350//
351//
353{
354 ResetTimer(0);
355
356 //
357 // Convert matrix to the form that Klu expects (Ap, Ai, VecAval)
358 //
359
360 if (MyPID_==0) {
361 assert( NumGlobalElements_ == SerialMatrix_->NumGlobalRows());
362 assert( NumGlobalElements_ == SerialMatrix_->NumGlobalCols());
363
364 if ( ! AddZeroToDiag_ ) {
365 assert( numentries_ == SerialMatrix_->NumGlobalNonzeros()) ;
366 } else {
367 numentries_ = SerialMatrix_->NumGlobalNonzeros() ;
368 }
369
370 Epetra_CrsMatrix *CrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(SerialMatrix_);
371 bool StorageOptimized = ( CrsMatrix != 0 && CrsMatrix->StorageOptimized() );
372
373 if ( AddToDiag_ != 0.0 ) StorageOptimized = false ;
374
375 if ( firsttime ) {
376 Ap.resize( NumGlobalElements_+1 );
377 Ai.resize( EPETRA_MAX( NumGlobalElements_, numentries_) ) ;
378 if ( ! StorageOptimized ) {
379 VecAval.resize( EPETRA_MAX( NumGlobalElements_, numentries_) ) ;
380 Aval = &VecAval[0];
381 }
382 }
383
384 double *RowValues;
385 int *ColIndices;
386 int NumEntriesThisRow;
387
388 if( StorageOptimized ) {
389 if ( firsttime ) {
390 Ap[0] = 0;
391 for ( int MyRow = 0; MyRow <NumGlobalElements_; MyRow++ ) {
392 if( CrsMatrix->
393 ExtractMyRowView( MyRow, NumEntriesThisRow, RowValues,
394 ColIndices ) != 0 )
395 AMESOS_CHK_ERR( -6 );
396 for ( int j=0; j < NumEntriesThisRow; j++ ) {
397 Ai[Ap[MyRow]+j] = ColIndices[j];
398 }
399 if ( MyRow == 0 ) {
400 Aval = RowValues ;
401 }
402 Ap[MyRow+1] = Ap[MyRow] + NumEntriesThisRow ;
403 }
404 }
405 } else {
406
407 int Ai_index = 0 ;
408 int MyRow;
409
410 int MaxNumEntries_ = SerialMatrix_->MaxNumEntries();
411 if ( firsttime && CrsMatrix == 0 ) {
412 ColIndicesV_.resize(MaxNumEntries_);
413 RowValuesV_.resize(MaxNumEntries_);
414 }
415
416 for ( MyRow = 0; MyRow <NumGlobalElements_; MyRow++ ) {
417 if ( CrsMatrix != 0 ) {
418 if( CrsMatrix->
419 ExtractMyRowView( MyRow, NumEntriesThisRow, RowValues,
420 ColIndices ) != 0 )
421 AMESOS_CHK_ERR( -6 );
422
423 } else {
424 if( SerialMatrix_->
425 ExtractMyRowCopy( MyRow, MaxNumEntries_,
426 NumEntriesThisRow, &RowValuesV_[0],
427 &ColIndicesV_[0] ) != 0 )
428 AMESOS_CHK_ERR( -6 );
429
430 RowValues = &RowValuesV_[0];
431 ColIndices = &ColIndicesV_[0];
432 }
433
434 if ( firsttime ) {
435 Ap[MyRow] = Ai_index ;
436 for ( int j = 0; j < NumEntriesThisRow; j++ ) {
437 Ai[Ai_index] = ColIndices[j] ;
438 // VecAval[Ai_index] = RowValues[j] ; // We have to do this because of the hacks to get around bug #1502
439 if (ColIndices[j] == MyRow) {
440 VecAval[Ai_index] += AddToDiag_;
441 }
442 Ai_index++;
443 }
444 } else {
445 for ( int j = 0; j < NumEntriesThisRow; j++ ) {
446 VecAval[Ai_index] = RowValues[j] ;
447 if (ColIndices[j] == MyRow) {
448 VecAval[Ai_index] += AddToDiag_;
449 }
450 Ai_index++;
451 }
452 }
453 }
454 Ap[MyRow] = Ai_index ;
455 }
457 cholmod_sparse& pk_A = PrivateParakleteData_->pk_A_ ;
458 pk_A.nrow = NumGlobalElements_ ;
459 pk_A.ncol = NumGlobalElements_ ;
460 pk_A.nzmax = numentries_ ;
461 pk_A.p = &Ap[0] ;
462 pk_A.i = &Ai[0] ;
463 pk_A.nz = 0;
464 if ( firsttime ) {
465 pk_A.x = 0 ;
466 pk_A.xtype = CHOLMOD_PATTERN ;
467 }
468 else
469 {
470 pk_A.x = Aval ;
471 pk_A.xtype = CHOLMOD_REAL ;
472 }
473
474 pk_A.z = 0 ;
475 pk_A.stype = 0 ; // symmetric
476 pk_A.itype = CHOLMOD_LONG ;
477 pk_A.dtype = CHOLMOD_DOUBLE ;
478 pk_A.sorted = 0 ;
479 pk_A.packed = 1 ;
480 }
481
482 MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
483
484 return 0;
485}
486
487//=============================================================================
488int Amesos_Paraklete::SetParameters( Teuchos::ParameterList &ParameterList ) {
489
490 // ========================================= //
491 // retrive PARAKLETE's parameters from list. //
492 // default values defined in the constructor //
493 // ========================================= //
494
495 // retrive general parameters
496
497 SetStatusParameters( ParameterList );
498 SetControlParameters( ParameterList );
499
500
501
502 if (ParameterList.isParameter("TrustMe") )
503 TrustMe_ = ParameterList.get<bool>( "TrustMe" );
504
505#if 0
506
507 unused for now
508
509 if (ParameterList.isSublist("Paraklete") ) {
510 Teuchos::ParameterList ParakleteParams = ParameterList.sublist("Paraklete") ;
511 }
512#endif
513
514 return 0;
515}
516
517
518//=============================================================================
520{
521 ResetTimer(0);
522
523 if ( IamInGroup_ )
524 PrivateParakleteData_->LUsymbolic_ =
525 rcp( amesos_paraklete_analyze ( &PrivateParakleteData_->pk_A_, &*PrivateParakleteData_->common_ )
526 ,deallocFunctorDeleteWithCommon<paraklete_symbolic>(PrivateParakleteData_->common_,
527 amesos_paraklete_free_symbolic)
528 ,true
529 );
530
531 SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_, 0);
532
533 return 0;
534}
535
536//=============================================================================
538{
539 // Changed this; it was "if (!TrustMe)...
540 // The behavior is not intuitive. Maybe we should introduce a new
541 // parameter, FastSolvers or something like that, that does not perform
542 // any AddTime, ResetTimer, GetTime.
543
544 ResetTimer(0);
545
546 //bool factor_with_pivoting = true ;
547
548 if ( IamInGroup_ )
549 PrivateParakleteData_->LUnumeric_ =
550 rcp( amesos_paraklete_factorize ( &PrivateParakleteData_->pk_A_,
551 &*PrivateParakleteData_->LUsymbolic_,
552 &*PrivateParakleteData_->common_ )
553 ,deallocFunctorDeleteWithCommon<paraklete_numeric>(PrivateParakleteData_->common_,amesos_paraklete_free_numeric)
554 ,true
555 );
556
557 NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
558
559 return 0;
560}
561
562//=============================================================================
564 bool OK = true;
565
566 // Comment by Tim: The following code seems suspect. The variable "OK"
567 // is not set if the condition is true.
568 // Does the variable "OK" default to true?
569 if ( GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
570 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) {
571 OK = false;
572 }
573 return OK;
574}
575
576
577extern "C" {
578 void my_handler (int status, char *file, int line, char *msg)
579 {
580 printf ("Error handler: %s %d %d: %s\n", file, line, status, msg) ;
581 if (status != CHOLMOD_OK)
582 {
583 fprintf (stderr, "\n\n*********************************************\n");
584 fprintf (stderr, "**** Test failure: %s %d %d %s\n", file, line,
585 status, msg) ;
586 fprintf (stderr, "*********************************************\n\n");
587 fflush (stderr) ;
588 fflush (stdout) ;
589 }
590 }
591}
592
593//=============================================================================
595{
596 MyPID_ = Comm().MyPID();
597 NumProcs_ = Comm().NumProc();
598
601
602#ifdef HAVE_AMESOS_EPETRAEXT
603 transposer_ = static_cast<Teuchos::ENull>( 0 );
604#endif
605
606 CreateTimer(Comm(), 2);
607
608 ResetTimer(1);
609
610 // "overhead" time for the following method is considered here
612 assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols() );
613
614
616
617 //
618 // Perform checks in SymbolicFactorization(), but none in
619 // NumericFactorization() or Solve()
620 //
621 assert( ! TrustMe_ ) ;
622 if ( TrustMe_ ) {
623 if ( CrsMatrixA_ == 0 ) AMESOS_CHK_ERR(10 );
624 if( UseDataInPlace_ != 1 ) AMESOS_CHK_ERR( 10 ) ;
625 if( Reindex_ ) AMESOS_CHK_ERR( 10 ) ;
626 if( ! Problem_->GetLHS() ) AMESOS_CHK_ERR( 10 ) ;
627 if( ! Problem_->GetRHS() ) AMESOS_CHK_ERR( 10 ) ;
628 if( ! Problem_->GetLHS()->NumVectors() ) AMESOS_CHK_ERR( 10 ) ;
629 if( ! Problem_->GetRHS()->NumVectors() ) AMESOS_CHK_ERR( 10 ) ;
630 SerialB_ = Problem_->GetRHS() ;
631 SerialX_ = Problem_->GetLHS() ;
632 NumVectors_ = SerialX_->NumVectors();
633 if (MyPID_ == 0) {
636 }
637 }
638
639
640 PrivateParakleteData_->common_ = rcp(new paraklete_common());
641
642 const Epetra_MpiComm* MpiComm = dynamic_cast<const Epetra_MpiComm*>(&Comm());
643 assert (MpiComm != 0);
644
645 MPI_Comm PK_Comm;
646 //
647 // Create an MPI group with MaxProcesses_ processes
648 //
649 if ( MaxProcesses_ != Comm().NumProc()) {
650 if(ParakleteComm_) {
651 MPI_Comm_free(&ParakleteComm_);
652 ParakleteComm_ = 0 ;
653 }
654 std::vector<int> ProcsInGroup(MaxProcesses_);
655 IamInGroup_ = false;
656 for (int i = 0 ; i < MaxProcesses_ ; ++i) {
657 ProcsInGroup[i] = i;
658 if ( Comm().MyPID() == i ) IamInGroup_ = true;
659 }
660
661 MPI_Group OrigGroup, ParakleteGroup;
662 MPI_Comm_group(MpiComm->GetMpiComm(), &OrigGroup);
663 MPI_Group_incl(OrigGroup, MaxProcesses_, &ProcsInGroup[0], &ParakleteGroup);
664 MPI_Comm_create(MpiComm->GetMpiComm(), ParakleteGroup, &ParakleteComm_);
665 PK_Comm = ParakleteComm_ ;
666 } else {
667 IamInGroup_ = true;
668 PK_Comm = MpiComm->GetMpiComm() ;
669 }
670
671 paraklete_common& pk_common = *PrivateParakleteData_->common_ ;
672 cholmod_common *cm = &(pk_common.cm) ;
673 amesos_cholmod_l_start (cm) ;
674 PK_DEBUG_INIT ("pk", cm) ;
675 pk_common.nproc = MaxProcesses_ ;
676 pk_common.myid = Comm().MyPID() ;
677 //pk_common.mpicomm = PK_Comm ;
678 cm->print = 1 ;
679 cm->precise = TRUE ;
680 cm->error_handler = my_handler ;
681
682 pk_common.tol_diag = 0.001 ;
683 pk_common.tol_offdiag = 0.1 ;
684 pk_common.growth = 2. ;
685
686
687
688 // "overhead" time for the following two methods is considered here
690
692
693 OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
694
695 // All this time if PARAKLETE time
697
699
701
702 return 0;
703}
704
705//=============================================================================
707{
708 if ( true || !TrustMe_ )
709 {
711 if (IsSymbolicFactorizationOK_ == false)
713
714 ResetTimer(1); // "overhead" time
715
716 Epetra_CrsMatrix *CrsMatrixA = dynamic_cast<Epetra_CrsMatrix *>(RowMatrixA_);
717 if ( false && CrsMatrixA == 0 ) // hack to get around Bug #1502
719 assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols() );
720 if ( AddZeroToDiag_ == 0 )
721 assert( numentries_ == RowMatrixA_->NumGlobalNonzeros() );
722
724
725 if ( false && CrsMatrixA == 0 ) { // continuation of hack to avoid bug #1502
727 } else {
729 }
730
731 OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
732 }
733
734 // this time is all for PARAKLETE
736
738
740
741 return 0;
742}
743
744//=============================================================================
746{
747 Epetra_MultiVector* vecX=0;
748 Epetra_MultiVector* vecB=0;
749
750 if ( !TrustMe_ ) {
751
752 SerialB_ = Problem_->GetRHS() ;
753 SerialX_ = Problem_->GetLHS() ;
754
755 Epetra_MultiVector* OrigVecX ;
756 Epetra_MultiVector* OrigVecB ;
757
758 if (IsNumericFactorizationOK_ == false)
760
761 ResetTimer(1);
762
763 //
764 // Reindex the LHS and RHS
765 //
766 OrigVecX = Problem_->GetLHS() ;
767 OrigVecB = Problem_->GetRHS() ;
768
769 if ( Reindex_ ) {
770#ifdef HAVE_AMESOS_EPETRAEXT
771 vecX = StdIndexDomain_->StandardizeIndex( OrigVecX ) ;
772 vecB = StdIndexRange_->StandardizeIndex( OrigVecB ) ;
773#else
774 AMESOS_CHK_ERR( -13 ) ; // Amesos_Paraklete can't handle non-standard indexing without EpetraExt
775#endif
776 } else {
777 vecX = OrigVecX ;
778 vecB = OrigVecB ;
779 }
780
781 if ((vecX == 0) || (vecB == 0))
782 AMESOS_CHK_ERR(-1); // something wrong in input
783
784 // Extract Serial versions of X and B
785
786 ResetTimer(0);
787
788 // Copy B to the serial version of B
789 //
790 if (false && UseDataInPlace_ == 1) {
791 SerialB_ = vecB;
792 SerialX_ = vecX;
793 NumVectors_ = Problem_->GetRHS()->NumVectors() ;
794 } else {
795 assert (UseDataInPlace_ == 0);
796
797 if( NumVectors_ != Problem_->GetRHS()->NumVectors() ) {
798 NumVectors_ = Problem_->GetRHS()->NumVectors() ;
799 SerialXextract_ = rcp( new Epetra_MultiVector(*SerialMap_,NumVectors_));
800 SerialBextract_ = rcp (new Epetra_MultiVector(*SerialMap_,NumVectors_));
801 }
802 if (NumVectors_ != vecB->NumVectors())
803 AMESOS_CHK_ERR(-1); // internal error
804
805 ImportRangeToSerial_ = rcp(new Epetra_Import ( *SerialMap_, vecB->Map() ) );
806 if ( SerialBextract_->Import(*vecB,*ImportRangeToSerial_,Insert) )
807 AMESOS_CHK_ERR( -1 ) ; // internal error
808
811 }
812 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
813
814 // Call PARAKLETE to perform the solve
815
816 ResetTimer(0);
817 if (MyPID_ == 0) {
821 AMESOS_CHK_ERR(-1);
822 }
823
824 OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
825 }
826 if ( MyPID_ == 0) {
827 if ( NumVectors_ == 1 ) {
828 for ( int i = 0 ; i < NumGlobalElements_ ; i++ )
830 } else {
831 SerialX_->Scale(1.0, *SerialB_ ) ; // X = B (Klu overwrites B with X)
832 }
833 }
834 if ( IamInGroup_ )
835 for (int i = 0; i < NumVectors_ ; i++ ) {
836 amesos_paraklete_solve( &*PrivateParakleteData_->LUnumeric_, &*PrivateParakleteData_->LUsymbolic_,
838 }
839 SolveTime_ = AddTime("Total solve time", SolveTime_, 0);
840
841 // Copy X back to the original vector
842
843 ResetTimer(0);
844 ResetTimer(1);
845
846 if (UseDataInPlace_ == 0) {
847 ImportDomainToSerial_ = rcp(new Epetra_Import ( *SerialMap_, vecX->Map() ) );
848 vecX->Export( *SerialX_, *ImportDomainToSerial_, Insert ) ;
849
850 } // otherwise we are already in place.
851
852 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
853
854#if 0
855 //
856 // ComputeTrueResidual causes TestOptions to fail on my linux box // Bug #1417
858 ComputeTrueResidual(*SerialMatrix_, *vecX, *vecB, UseTranspose(), "Amesos_Paraklete");
859#endif
860
862 ComputeVectorNorms(*vecX, *vecB, "Amesos_Paraklete");
863
864 OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
865
866 ++NumSolve_;
867
868 return(0) ;
869}
870
871// ================================================ ====== ==== ==== == =
872
874{
875
876 if (MyPID_) return;
877
878 PrintLine();
879
880 std::cout << "Amesos_Paraklete : Matrix has " << NumGlobalElements_ << " rows"
881 << " and " << numentries_ << " nonzeros" << std::endl;
882 std::cout << "Amesos_Paraklete : Nonzero elements per row = "
883 << 1.0*numentries_/NumGlobalElements_ << std::endl;
884 std::cout << "Amesos_Paraklete : Percentage of nonzero elements = "
885 << 100.0*numentries_/(pow(double(NumGlobalElements_),double(2.0))) << std::endl;
886 std::cout << "Amesos_Paraklete : Use transpose = " << UseTranspose_ << std::endl;
887
888 PrintLine();
889
890 return;
891
892}
893
894// ================================================ ====== ==== ==== == =
895
897{
898 if (MyPID_) return;
899
900 double ConTime = GetTime(MtxConvTime_);
901 double MatTime = GetTime(MtxRedistTime_);
902 double VecTime = GetTime(VecRedistTime_);
903 double SymTime = GetTime(SymFactTime_);
904 double NumTime = GetTime(NumFactTime_);
905 double SolTime = GetTime(SolveTime_);
906 double OveTime = GetTime(OverheadTime_);
907
909 SymTime /= NumSymbolicFact_;
910
911 if (NumNumericFact_)
912 NumTime /= NumNumericFact_;
913
914 if (NumSolve_)
915 SolTime /= NumSolve_;
916
917 std::string p = "Amesos_Paraklete : ";
918 PrintLine();
919
920 std::cout << p << "Time to convert matrix to Paraklete format = "
921 << ConTime << " (s)" << std::endl;
922 std::cout << p << "Time to redistribute matrix = "
923 << MatTime << " (s)" << std::endl;
924 std::cout << p << "Time to redistribute vectors = "
925 << VecTime << " (s)" << std::endl;
926 std::cout << p << "Number of symbolic factorizations = "
927 << NumSymbolicFact_ << std::endl;
928 std::cout << p << "Time for sym fact = "
929 << SymTime * NumSymbolicFact_ << " (s), avg = " << SymTime << " (s)" << std::endl;
930 std::cout << p << "Number of numeric factorizations = "
931 << NumNumericFact_ << std::endl;
932 std::cout << p << "Time for num fact = "
933 << NumTime * NumNumericFact_ << " (s), avg = " << NumTime << " (s)" << std::endl;
934 std::cout << p << "Number of solve phases = "
935 << NumSolve_ << std::endl;
936 std::cout << p << "Time for solve = "
937 << SolTime * NumSolve_ << " (s), avg = " << SolTime << " (s)" << std::endl;
938
939 double tt = SymTime * NumSymbolicFact_ + NumTime * NumNumericFact_ + SolTime * NumSolve_;
940 if (tt != 0)
941 {
942 std::cout << p << "Total time spent in Amesos = " << tt << " (s) " << std::endl;
943 std::cout << p << "Total time spent in the Amesos interface = " << OveTime << " (s)" << std::endl;
944 std::cout << p << "(the above time does not include PARAKLETE time)" << std::endl;
945 std::cout << p << "Amesos interface time / total time = " << OveTime / tt << std::endl;
946 }
947
948 PrintLine();
949
950 return;
951}
952
953
#define AMESOS_CHK_ERR(a)
void my_handler(int status, char *file, int line, char *msg)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
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< paraklete_symbolic > LUsymbolic_
Teuchos::RCP< paraklete_common > common_
Teuchos::RCP< paraklete_numeric > LUnumeric_
bool UseTranspose() const
Returns the current UseTranspose setting.
bool MatrixShapeOK() const
Returns true if PARAKLETE can handle this matrix shape.
bool TrustMe_
If true, no checks are made and the matrix is assume to be distributed.
Epetra_MultiVector * SerialX_
std::vector< int > ColIndicesV_
Only used for RowMatrices to extract copies.
std::vector< long > Ai
int MtxConvTime_
Quick access pointers to internal timing information.
int UseDataInPlace_
1 if Problem_->GetOperator() is stored entirely on process 0
Teuchos::RCP< Epetra_MultiVector > SerialXextract_
Serial versions of the LHS and RHS (if necessary)
int numentries_
Number of non-zero entries in Problem_->GetOperator()
std::vector< long > Ap
Ap, Ai, Aval form the compressed row storage used by Paraklete Ai and Aval can point directly into a ...
Epetra_MultiVector * SerialB_
Serial versions of the LHS and RHS (may point to the original vector if serial)
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrixA_
Points to a Serial Copy of A (unused if UseDataInPlace_==1)
Epetra_RowMatrix * SerialMatrix_
Points to a Serial Copy of A.
void PrintStatus() const
Prints information about the factorization and solution phases.
const Epetra_LinearProblem * GetProblem() const
Get a pointer to the Problem.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
std::vector< double > RowValuesV_
Only used for RowMatrices to extract copies.
Teuchos::RCP< Amesos_StandardIndex > StdIndexDomain_
int NumVectors_
Number of vectors in RHS and LHS.
double * SerialXBvalues_
Pointer to the actual values in the serial version of X and B.
Teuchos::RCP< Amesos_StandardIndex > StdIndex_
void PrintTiming() const
Prints timing information.
Epetra_CrsMatrix * CrsMatrixA_
Operator converted to a CrsMatrix.
Teuchos::RCP< Amesos_Paraklete_Pimpl > PrivateParakleteData_
Teuchos::RCP< Amesos_StandardIndex > StdIndexRange_
int ConvertToParakleteCRS(bool firsttime)
Teuchos::RCP< Epetra_MultiVector > SerialBextract_
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
int Solve()
Solves A X = B (or AT x = B)
const Epetra_LinearProblem * Problem_
Pointer to the linear system problem.
bool UseTranspose_
If true, the transpose of A is used.
Epetra_RowMatrix * StdIndexMatrix_
Points to a Contiguous Copy of A.
~Amesos_Paraklete(void)
Amesos_Paraklete Destructor.
Teuchos::RCP< Epetra_Map > SerialMap_
Points to a Serial Map (unused if UseDataInPlace_ == 1 )
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
Amesos_Paraklete(const Epetra_LinearProblem &LinearProblem)
Amesos_Paraklete Constructor.
std::vector< double > VecAval
Teuchos::RCP< Epetra_Import > ImportDomainToSerial_
int NumGlobalElements_
Number of rows and columns in the Problem_->GetOperator()
Teuchos::RCP< Epetra_Import > ImportRangeToSerial_
Epetra_RowMatrix * RowMatrixA_
Operator converted to a RowMatrix.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Teuchos::RCP< Epetra_Import > ImportToSerial_
Importer to process 0.
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 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.