33 #include "Epetra_Map.h" 34 #include "Epetra_Import.h" 35 #include "Epetra_Export.h" 36 #include "Epetra_RowMatrix.h" 37 #include "Epetra_CrsMatrix.h" 38 #include "Epetra_Vector.h" 39 #include "Epetra_Util.h" 57 RcondValidOnAllProcs_(true),
80 amd_defaults(control);
106 const Epetra_Map &OriginalMap =
Matrix()->RowMatrixRowMap() ;
112 int NumMyElements_ = 0 ;
115 IsLocal_ = ( OriginalMap.NumMyElements() ==
116 OriginalMap.NumGlobalElements() )?1:0;
152 I was not able to make
this work - 11 Feb 2006
156 SerialCrsMatrix().SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) ) ;
162 for (
int i = 0 ; i <
SerialMap_->NumGlobalElements(); i++ )
199 int NumEntriesThisRow;
205 Ap[MyRow] = Ai_index ;
208 &
Aval[Ai_index], &
Ai[Ai_index]);
215 for (
int i = 0 ; i < NumEntriesThisRow ; ++i) {
216 if (
Ai[Ai_index+i] == MyRow) {
223 Ai_index += NumEntriesThisRow;
226 Ap[MyRow] = Ai_index ;
258 double *Control = (
double *) NULL, *Info = (
double *) NULL;
261 umfpack_di_free_symbolic (&
Symbolic) ;
266 symbolic_ok = (status == UMFPACK_OK);
270 Comm().Broadcast(&symbolic_ok, 1, 0);
295 std::vector<double> Control(UMFPACK_CONTROL);
296 std::vector<double> Info(UMFPACK_INFO);
297 umfpack_di_defaults( &Control[0] ) ;
299 int status = umfpack_di_numeric (&
Ap[0],
306 Rcond_ = Info[UMFPACK_RCOND];
309 std::cout <<
" Rcond_ = " <<
Rcond_ << std::endl ;
314 int * Lp = (
int *) malloc ((
n+1) *
sizeof (int)) ;
315 int * Lj = (
int *) malloc (lnz1 *
sizeof (
int)) ;
316 double * Lx = (
double *) malloc (lnz1 *
sizeof (
double)) ;
317 int * Up = (
int *) malloc ((
n+1) *
sizeof (int)) ;
318 int * Ui = (
int *) malloc (unz1 *
sizeof (
int)) ;
319 double * Ux = (
double *) malloc (unz1 *
sizeof (
double)) ;
320 int * P = (
int *) malloc (
n *
sizeof (
int)) ;
321 int * Q = (
int *) malloc (
n *
sizeof (
int)) ;
322 double * Dx = (
double *) NULL ;
323 double * Rs = (
double *) malloc (
n *
sizeof (
double)) ;
324 if (!Lp || !Lj || !Lx || !Up || !Ui || !Ux || !P || !Q || !Rs)
329 status = umfpack_di_get_numeric (Lp, Lj, Lx, Up, Ui, Ux,
330 P, Q, Dx, &do_recip, Rs,
Numeric) ;
336 printf (
"\nL (lower triangular factor of C): ") ;
337 (void) umfpack_di_report_matrix (
n,
n, Lp, Lj, Lx, 0, &Control[0]) ;
338 printf (
"\nU (upper triangular factor of C): ") ;
339 (void) umfpack_di_report_matrix (
n,
n, Up, Ui, Ux, 1, &Control[0]) ;
341 (void) umfpack_di_report_perm (
n, P, &Control[0]) ;
343 (void) umfpack_di_report_perm (
n, Q, &Control[0]) ;
344 printf (
"\nScale factors: row i of A is to be ") ;
348 numeric_ok = (status == UMFPACK_OK);
352 Comm().Broadcast(&numeric_ok, 1, 0);
382 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
383 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) OK =
false;
478 Epetra_MultiVector* vecX =
Problem_->GetLHS();
479 Epetra_MultiVector* vecB =
Problem_->GetRHS();
481 if ((vecX == 0) || (vecB == 0))
484 int NumVectors = vecX->NumVectors();
485 if (NumVectors != vecB->NumVectors())
488 Epetra_MultiVector *SerialB, *SerialX;
492 double *SerialXvalues ;
493 double *SerialBvalues ;
495 Epetra_MultiVector* SerialXextract = 0;
496 Epetra_MultiVector* SerialBextract = 0;
507 SerialXextract =
new Epetra_MultiVector(
SerialMap(),NumVectors);
508 SerialBextract =
new Epetra_MultiVector(
SerialMap(),NumVectors);
510 SerialBextract->Import(*vecB,
Importer(),Insert);
511 SerialB = SerialBextract;
512 SerialX = SerialXextract;
525 int SerialBlda, SerialXlda ;
526 int UmfpackRequest =
UseTranspose()?UMFPACK_A:UMFPACK_At ;
531 ierr = SerialB->ExtractView(&SerialBvalues, &SerialBlda);
533 ierr = SerialX->ExtractView(&SerialXvalues, &SerialXlda);
538 for (
int j =0 ; j < NumVectors; j++ ) {
539 double *Control = (
double *) NULL, *Info = (
double *) NULL ;
541 status = umfpack_di_solve (UmfpackRequest, &
Ap[0],
543 &SerialXvalues[j*SerialXlda],
544 &SerialBvalues[j*SerialBlda],
559 vecX->Export(*SerialX,
Importer(), Insert ) ;
560 if (SerialBextract)
delete SerialBextract ;
561 if (SerialXextract)
delete SerialXextract ;
568 Epetra_RowMatrix*
Matrix =
569 dynamic_cast<Epetra_RowMatrix*
>(
Problem_->GetOperator());
592 <<
" and " <<
numentries_ <<
" nonzeros" << std::endl;
593 std::cout <<
"Amesos_Umfpack : Nonzero elements per row = " 595 std::cout <<
"Amesos_Umfpack : Percentage of nonzero elements = " 597 std::cout <<
"Amesos_Umfpack : Use transpose = " <<
UseTranspose_ << std::endl;
627 std::string p =
"Amesos_Umfpack : ";
630 std::cout << p <<
"Time to convert matrix to Umfpack format = " 631 << ConTime <<
" (s)" << std::endl;
632 std::cout << p <<
"Time to redistribute matrix = " 633 << MatTime <<
" (s)" << std::endl;
634 std::cout << p <<
"Time to redistribute vectors = " 635 << VecTime <<
" (s)" << std::endl;
636 std::cout << p <<
"Number of symbolic factorizations = " 638 std::cout << p <<
"Time for sym fact = " 640 << SymTime <<
" (s)" << std::endl;
641 std::cout << p <<
"Number of numeric factorizations = " 643 std::cout << p <<
"Time for num fact = " 645 << NumTime <<
" (s)" << std::endl;
646 std::cout << p <<
"Number of solve phases = " 648 std::cout << p <<
"Time for solve = " 649 << SolTime *
NumSolve_ <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
653 std::cout << p <<
"Total time spent in Amesos = " << tt <<
" (s) " << std::endl;
654 std::cout << p <<
"Total time spent in the Amesos interface = " << OveTime <<
" (s)" << std::endl;
655 std::cout << p <<
"(the above time does not include UMFPACK time)" << std::endl;
656 std::cout << p <<
"Amesos interface time / total time = " << OveTime / tt << std::endl;
std::vector< int > Ap
Ap, Ai, Aval form the compressed row storage used by Umfpack.
int numentries_
Number of non-zero entries in Problem_->GetOperator()
int NumSymbolicFact_
Number of symbolic factorization phases.
Amesos_Klu: A serial, unblocked code ideal for getting started and for very sparse matrices...
double GetTime(const std::string what) const
Gets the cumulative time using the string.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
int PerformSymbolicFactorization()
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
int MtxConvTime_
Quick access pointers to internal timer data.
const Epetra_Map & SerialMap() const
const Epetra_LinearProblem * Problem_
Pointer to the linear problem to solve.
void * Numeric
Umfpack internal opaque object.
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
bool MatrixShapeOK() const
Returns true if UMFPACK can handle this matrix shape.
const Epetra_CrsMatrix & SerialCrsMatrix() const
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrixA_
bool AddZeroToDiag_
Adds zero to diagonal of redistributed matrix (some solvers choke on a matrix with a partly empty dia...
void PrintStatus() const
Prints information about the factorization and solution phases.
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
double GetRcond() const
Returns an estimate of the reciprocal of the condition number.
int NumNumericFact_
Number of numeric factorization phases.
~Amesos_Umfpack(void)
Amesos_Umfpack Destructor.
int ConvertToSerial(const bool FirstTime)
Converts matrix to a serial Epetra_CrsMatrix.
int NumSolve_
Number of solves.
int NumGlobalElements_
Number of rows and columns in the Problem_->GetOperator()
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
#define AMESOS_CHK_ERR(a)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
Epetra_RowMatrix * Matrix()
Returns a pointer to the linear system matrix.
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.
bool PrintTiming_
If true, prints timing information in the destructor.
bool PrintStatus_
If true, print additional information in the destructor.
Teuchos::RCP< Epetra_Import > ImportToSerial_
Importer from distributed to serial (all rows on process 0).
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
Epetra_RowMatrix * SerialMatrix_
Points to a Serial Copy of A.
std::vector< double > Aval
Amesos_Umfpack(const Epetra_LinearProblem &LinearProblem)
Amesos_Umfpack Constructor.
const Epetra_Import & Importer() const
int ConvertToUmfpackCRS()
bool UseTranspose() const
Returns the current UseTranspose setting.
double Rcond_
Reciprocal condition number estimate.
int Solve()
Solves A X = B (or AT x = B)
void PrintTiming() const
Prints timing information.
int PerformNumericFactorization()
const int NumericallySingularMatrixError
int NumericFactorization()
Performs NumericFactorization on the matrix A.
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.
const int StructurallySingularMatrixError
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
bool UseTranspose_
If true, solve the problem with the transpose.
int IsLocal_
1 if Problem_->GetOperator() is stored entirely on process 0
int verbose_
Toggles the output level.
bool RcondValidOnAllProcs_
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
void PrintLine() const
Prints line on std::cout.
Teuchos::RCP< Epetra_Map > SerialMap_
Points to a Serial Map (unused if IsLocal == 1 )
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
void * Symbolic
Umfpack internal opaque object.
double AddToDiag_
Add this value to the diagonal.