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 ;
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;
const int NumericallySingularMatrixError
const int StructurallySingularMatrixError
#define AMESOS_CHK_ERR(a)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
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.
Amesos_Klu: A serial, unblocked code ideal for getting started and for very sparse matrices,...
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.
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
bool UseTranspose_
If true, solve the problem with the transpose.
int IsLocal_
1 if Problem_->GetOperator() is stored entirely on process 0
bool UseTranspose() const
Returns the current UseTranspose setting.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Teuchos::RCP< Epetra_Import > ImportToSerial_
Importer from distributed to serial (all rows on process 0).
int PerformNumericFactorization()
void PrintTiming() const
Prints timing information.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
bool MatrixShapeOK() const
Returns true if UMFPACK can handle this matrix shape.
const Epetra_LinearProblem * Problem_
Pointer to the linear problem to solve.
void * Numeric
Umfpack internal opaque object.
Amesos_Umfpack(const Epetra_LinearProblem &LinearProblem)
Amesos_Umfpack Constructor.
int PerformSymbolicFactorization()
const Epetra_CrsMatrix & SerialCrsMatrix() const
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrixA_
std::vector< double > Aval
bool RcondValidOnAllProcs_
int ConvertToUmfpackCRS()
const Epetra_Map & SerialMap() const
void PrintStatus() const
Prints information about the factorization and solution phases.
const Epetra_Import & Importer() const
~Amesos_Umfpack(void)
Amesos_Umfpack Destructor.
int numentries_
Number of non-zero entries in Problem_->GetOperator()
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
void * Symbolic
Umfpack internal opaque object.
std::vector< int > Ap
Ap, Ai, Aval form the compressed row storage used by Umfpack.
int Solve()
Solves A X = B (or AT x = B)
Epetra_RowMatrix * Matrix()
Returns a pointer to the linear system matrix.
int NumGlobalElements_
Number of rows and columns in the Problem_->GetOperator()
Teuchos::RCP< Epetra_Map > SerialMap_
Points to a Serial Map (unused if IsLocal == 1 )
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
int ConvertToSerial(const bool FirstTime)
Converts matrix to a serial Epetra_CrsMatrix.
double Rcond_
Reciprocal condition number estimate.
double GetRcond() const
Returns an estimate of the reciprocal of the condition number.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int MtxConvTime_
Quick access pointers to internal timer data.
Epetra_RowMatrix * SerialMatrix_
Points to a Serial Copy of 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.