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),
69 Teuchos::ParameterList ParamList ;
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 ;
256 double *Control = (
double *)
NULL, *Info = (
double *)
NULL;
259 umfpack_di_free_symbolic (&
Symbolic) ;
279 std::vector<double> Control(UMFPACK_CONTROL);
280 std::vector<double> Info(UMFPACK_INFO);
281 umfpack_di_defaults( &Control[0] ) ;
283 int status = umfpack_di_numeric (&
Ap[0],
290 Rcond_ = Info[UMFPACK_RCOND];
293 std::cout <<
" Rcond_ = " <<
Rcond_ << std::endl ;
298 int * Lp = (
int *) malloc ((n+1) *
sizeof (int)) ;
299 int * Lj = (
int *) malloc (lnz1 *
sizeof (
int)) ;
300 double * Lx = (
double *) malloc (lnz1 *
sizeof (
double)) ;
301 int * Up = (
int *) malloc ((n+1) *
sizeof (int)) ;
302 int * Ui = (
int *) malloc (unz1 *
sizeof (
int)) ;
303 double * Ux = (
double *) malloc (unz1 *
sizeof (
double)) ;
304 int *
P = (
int *) malloc (n *
sizeof (
int)) ;
305 int * Q = (
int *) malloc (n *
sizeof (
int)) ;
306 double * Dx = (
double *)
NULL ;
307 double * Rs = (
double *) malloc (n *
sizeof (
double)) ;
308 if (!Lp || !Lj || !Lx || !Up || !Ui || !Ux || !
P || !Q || !Rs)
313 status = umfpack_di_get_numeric (Lp, Lj, Lx, Up, Ui, Ux,
320 printf (
"\nL (lower triangular factor of C): ") ;
321 (void) umfpack_di_report_matrix (n, n, Lp, Lj, Lx, 0, &Control[0]) ;
322 printf (
"\nU (upper triangular factor of C): ") ;
323 (void) umfpack_di_report_matrix (n, n, Up, Ui, Ux, 1, &Control[0]) ;
325 (void) umfpack_di_report_perm (n,
P, &Control[0]) ;
327 (void) umfpack_di_report_perm (n, Q, &Control[0]) ;
328 printf (
"\nScale factors: row i of A is to be ") ;
332 assert( status == 0 ) ;
355 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
356 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) OK =
false;
451 Epetra_MultiVector* vecX =
Problem_->GetLHS();
452 Epetra_MultiVector* vecB =
Problem_->GetRHS();
454 if ((vecX == 0) || (vecB == 0))
457 int NumVectors = vecX->NumVectors();
458 if (NumVectors != vecB->NumVectors())
461 Epetra_MultiVector *SerialB, *SerialX;
465 double *SerialXvalues ;
466 double *SerialBvalues ;
468 Epetra_MultiVector* SerialXextract = 0;
469 Epetra_MultiVector* SerialBextract = 0;
480 SerialXextract =
new Epetra_MultiVector(
SerialMap(),NumVectors);
481 SerialBextract =
new Epetra_MultiVector(
SerialMap(),NumVectors);
483 SerialBextract->Import(*vecB,
Importer(),Insert);
484 SerialB = SerialBextract;
485 SerialX = SerialXextract;
498 int SerialBlda, SerialXlda ;
499 int UmfpackRequest =
UseTranspose()?UMFPACK_A:UMFPACK_At ;
504 ierr = SerialB->ExtractView(&SerialBvalues, &SerialBlda);
506 ierr = SerialX->ExtractView(&SerialXvalues, &SerialXlda);
511 for (
int j =0 ; j < NumVectors; j++ ) {
512 double *Control = (
double *)
NULL, *Info = (
double *)
NULL ;
514 status = umfpack_di_solve (UmfpackRequest, &
Ap[0],
516 &SerialXvalues[j*SerialXlda],
517 &SerialBvalues[j*SerialBlda],
532 vecX->Export(*SerialX,
Importer(), Insert ) ;
533 if (SerialBextract)
delete SerialBextract ;
534 if (SerialXextract)
delete SerialXextract ;
541 Epetra_RowMatrix*
Matrix =
542 dynamic_cast<Epetra_RowMatrix*
>(
Problem_->GetOperator());
565 <<
" and " <<
numentries_ <<
" nonzeros" << std::endl;
566 std::cout <<
"Amesos_Umfpack : Nonzero elements per row = " 568 std::cout <<
"Amesos_Umfpack : Percentage of nonzero elements = " 570 std::cout <<
"Amesos_Umfpack : Use transpose = " <<
UseTranspose_ << std::endl;
600 std::string p =
"Amesos_Umfpack : ";
603 std::cout << p <<
"Time to convert matrix to Umfpack format = " 604 << ConTime <<
" (s)" << std::endl;
605 std::cout << p <<
"Time to redistribute matrix = " 606 << MatTime <<
" (s)" << std::endl;
607 std::cout << p <<
"Time to redistribute vectors = " 608 << VecTime <<
" (s)" << std::endl;
609 std::cout << p <<
"Number of symbolic factorizations = " 611 std::cout << p <<
"Time for sym fact = " 613 << SymTime <<
" (s)" << std::endl;
614 std::cout << p <<
"Number of numeric factorizations = " 616 std::cout << p <<
"Time for num fact = " 618 << NumTime <<
" (s)" << std::endl;
619 std::cout << p <<
"Number of solve phases = " 621 std::cout << p <<
"Time for solve = " 622 << SolTime *
NumSolve_ <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
626 std::cout << p <<
"Total time spent in Amesos = " << tt <<
" (s) " << std::endl;
627 std::cout << p <<
"Total time spent in the Amesos interface = " << OveTime <<
" (s)" << std::endl;
628 std::cout << p <<
"(the above time does not include UMFPACK time)" << std::endl;
629 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)
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()
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.
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.