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" 116 Epetra_RowMatrix *RowMatrixA =
Problem_->GetMatrix();
120 const Epetra_Map& OriginalMap = RowMatrixA->RowMatrixRowMap() ;
121 const Epetra_MpiComm& comm1 =
dynamic_cast<const Epetra_MpiComm &
> (
Comm());
122 int numrows = RowMatrixA->NumGlobalRows();
123 int numentries = RowMatrixA->NumGlobalNonzeros();
127 Epetra_CrsMatrix* CastCrsMatrixA =
128 dynamic_cast<Epetra_CrsMatrix*
>(RowMatrixA);
132 Graph =
Teuchos::rcp(const_cast<Epetra_CrsGraph*>(&(CastCrsMatrixA->Graph())),
false);
136 int MaxNumEntries = RowMatrixA->MaxNumEntries();
137 Graph =
Teuchos::rcp(
new Epetra_CrsGraph(
Copy, OriginalMap, MaxNumEntries));
139 std::vector<int> Indices(MaxNumEntries);
140 std::vector<double> Values(MaxNumEntries);
142 for (
int i = 0 ; i < RowMatrixA->NumMyRows() ; ++i)
145 RowMatrixA->ExtractMyRowCopy(i, MaxNumEntries, NumEntries,
146 &Values[0], &Indices[0]);
148 for (
int j = 0 ; j < NumEntries ; ++j)
149 Indices[j] = RowMatrixA->RowMatrixColMap().GID(Indices[j]);
151 int GlobalRow = RowMatrixA->RowMatrixRowMap().GID(i);
152 Graph->InsertGlobalIndices(GlobalRow, NumEntries, &Indices[0]);
155 Graph->FillComplete();
161 std::vector<int> AllIDs( numrows ) ;
162 for (
int i = 0; i < numrows ; i++ ) AllIDs[i] = i ;
164 Epetra_Map ReplicatedMap( -1, numrows, &AllIDs[0], 0,
Comm());
165 Epetra_Import ReplicatedImporter(ReplicatedMap, OriginalMap);
166 Epetra_CrsGraph ReplicatedGraph(
Copy, ReplicatedMap, 0 );
168 AMESOS_CHK_ERR(ReplicatedGraph.Import(*Graph, ReplicatedImporter, Insert));
174 std::vector <int> Replicates(numrows);
175 std::vector <int> Ap(numrows + 1);
176 std::vector <int> Ai(EPETRA_MAX(numrows, numentries));
178 for(
int i = 0 ; i < numrows; i++ ) Replicates[i] = 1;
180 int NumEntriesPerRow ;
181 int *ColIndices = 0 ;
183 for (
int MyRow = 0; MyRow <numrows; MyRow++ ) {
184 AMESOS_CHK_ERR( ReplicatedGraph.ExtractMyRowView( MyRow, NumEntriesPerRow, ColIndices ) );
185 Ap[MyRow] = Ai_index ;
186 for (
int j = 0; j < NumEntriesPerRow; j++ ) {
187 Ai[Ai_index] = ColIndices[j] ;
191 assert( Ai_index == numentries ) ;
192 Ap[ numrows ] = Ai_index ;
202 std::vector<double> MyANonZ;
215 int NumGlobalNonzeros =
GetProblem()->GetMatrix()->NumGlobalNonzeros();
216 int NumRows =
GetProblem()->GetMatrix()->NumGlobalRows();
220 int OptNumProcs1 = 1+EPETRA_MAX( NumRows/10000, NumGlobalNonzeros/1000000 );
221 OptNumProcs1 = EPETRA_MIN(
NumProcs_,OptNumProcs1 );
225 int OptNumProcs2 = (int)sqrt(1.0 *
NumProcs_);
226 if( OptNumProcs2 < 1 ) OptNumProcs2 = 1;
258 int DscMax = DSC_Analyze( numrows, &Ap[0], &Ai[0], &Replicates[0] );
280 const int Limit = 5000000 ;
282 &MaxSingleBlock, Limit, DSC_LBLAS3, DSC_DBLAS2 ) ) ;
299 Epetra_RowMatrix* RowMatrixA =
Problem_->GetMatrix();
303 const Epetra_Map& OriginalMap = RowMatrixA->RowMatrixRowMap() ;
305 int numrows = RowMatrixA->NumGlobalRows();
306 assert( numrows == RowMatrixA->NumGlobalCols() );
311 std::vector<double> MyANonZ;
341 int max_num_entries = DscMat.MaxNumEntries() ;
342 std::vector<int> col_indices( max_num_entries ) ;
343 std::vector<double> mat_values( max_num_entries ) ;
350 typedef std::pair<int, double> Data;
351 std::vector<Data> sort_array(max_num_entries);
352 std::vector<int> sort_indices(max_num_entries);
356 int num_entries_this_row;
358 AMESOS_CHK_ERR( DscMat.ExtractMyRowCopy( i, max_num_entries, num_entries_this_row,
359 &mat_values[0], &col_indices[0] ) ) ;
361 AMESOS_CHK_ERR( DscMat.ExtractGlobalRowCopy( DscMat.GRID(i), max_num_entries, num_entries_this_row,
362 &mat_values[0], &col_indices[0] ) ) ;
374 for (
int j = 0; j < num_entries_this_row; j++ ) {
381 sort_array[j].second = mat_values[j] ;
383 sort(&sort_array[0], &sort_array[num_entries_this_row]);
385 for (
int j = 0; j < num_entries_this_row; j++ ) {
386 int NewColNumber = sort_array[j].first ;
387 if ( NewRowNumber <= NewColNumber ) MyANonZ[ NonZIndex++ ] = sort_array[j].second ;
398 const int SchemeCode = 1;
404 DSC_LLT, DSC_LBLAS3, DSC_DBLAS2 ) ) ;
417 bool OK =
GetProblem()->IsOperatorSymmetric() ;
423 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
424 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) OK =
false;
469 Epetra_RowMatrix *RowMatrixA =
Problem_->GetMatrix();
474 if (RowMatrixA->NumGlobalRows() != RowMatrixA->NumGlobalCols())
479 Epetra_MultiVector* vecX =
Problem_->GetLHS();
480 Epetra_MultiVector* vecB =
Problem_->GetRHS();
482 if ((vecX == 0) || (vecB == 0))
485 int NumVectors = vecX->NumVectors();
486 if (NumVectors != vecB->NumVectors())
489 double *dscmapXvalues ;
491 Epetra_MultiVector dscmapX(
DscRowMap(),NumVectors) ;
496 double *dscmapBvalues ;
498 Epetra_MultiVector dscmapB(
DscRowMap(), NumVectors ) ;
499 ierr = dscmapB.ExtractView( &dscmapBvalues, &dscmapBlda );
515 for (
int j =0 ; j < NumVectors; j++ ) {
532 vecX->Export( dscmapX,
Importer(), Insert ) ;
538 false,
"Amesos_Dscpack");
554 std::string p =
"Amesos_Dscpack : ";
557 int n =
GetProblem()->GetMatrix()->NumGlobalRows();
560 std::cout << p <<
"Matrix has " <<
n <<
" rows" 561 <<
" and " <<
nnz <<
" nonzeros" << std::endl;
562 std::cout << p <<
"Nonzero elements per row = " 563 << 1.0 *
nnz /
n << std::endl;
564 std::cout << p <<
"Percentage of nonzero elements = " 565 << 100.0 *
nnz /(pow(
n,2.0)) << std::endl;
566 std::cout << p <<
"Available process(es) = " <<
NumProcs_ << std::endl;
567 std::cout << p <<
"Process(es) used = " <<
DscNumProcs 569 std::cout << p <<
"Estimated total memory for factorization = " 605 std::string p =
"Amesos_Dscpack : ";
608 std::cout << p <<
"Time to convert matrix to DSCPACK format = " 609 << ConTime <<
" (s)" << std::endl;
610 std::cout << p <<
"Time to redistribute matrix = " 611 << MatTime <<
" (s)" << std::endl;
612 std::cout << p <<
"Time to redistribute vectors = " 613 << VecTime <<
" (s)" << std::endl;
614 std::cout << p <<
"Number of symbolic factorizations = " 616 std::cout << p <<
"Time for sym fact = " 617 << SymTime *
NumSymbolicFact_ <<
" (s), avg = " << SymTime <<
" (s)" << std::endl;
618 std::cout << p <<
"Number of numeric factorizations = " 620 std::cout << p <<
"Time for num fact = " 621 << NumTime *
NumNumericFact_ <<
" (s), avg = " << NumTime <<
" (s)" << std::endl;
622 std::cout << p <<
"Number of solve phases = " 624 std::cout << p <<
"Time for solve = " 625 << SolTime *
NumSolve_ <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
630 std::cout << p <<
"Total time spent in Amesos = " << tt <<
" (s) " << std::endl;
631 std::cout << p <<
"Total time spent in the Amesos interface = " << OveTime <<
" (s)" << std::endl;
632 std::cout << p <<
"(the above time does not include DSCPACK time)" << std::endl;
633 std::cout << p <<
"Amesos interface time / total time = " << OveTime / tt << std::endl;
const Epetra_Map & DscColMap() const
int NumSymbolicFact_
Number of symbolic factorization phases.
double GetTime(const std::string what) const
Gets the cumulative time using the string.
RCP< Epetra_Map > DscColMap_
RCP< Epetra_Import > Importer_
UF_long CHOLMOD() nnz(cholmod_sparse *A, cholmod_common *Common)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
bool A_and_LU_built
Tells us whether to free them.
void PrintStatus() const
Prints information about the factorization and solution phases.
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
const Epetra_Map & DscRowMap() const
int NumNumericFact_
Number of numeric factorization phases.
int CHOLMOD() sort(cholmod_sparse *A, cholmod_common *Common)
int NumSolve_
Number of solves.
const Epetra_LinearProblem * Problem_
Pointer to the linear problem.
Teuchos::RCP< Amesos_Dscpack_Pimpl > PrivateDscpackData_
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 SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
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.
int PerformSymbolicFactorization()
Performs the symbolic factorization.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
bool PrintStatus_
If true, print additional information in the destructor.
int * GlobalStructNewColNum
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
bool MatrixShapeOK() const
Returns true if DSCPACK can handle this matrix shape.
Amesos_Dscpack(const Epetra_LinearProblem &LinearProblem)
Amesos_Dscpack Constructor.
const Epetra_Import & Importer() const
~Amesos_Dscpack(void)
Amesos_Dscpack Destructor.
int Solve()
Solves A X = B (or AT x = B)
void PrintTiming() const
Prints timing information.
RCP< Epetra_Map > DscRowMap_
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.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
int PerformNumericFactorization()
Performs the numeric factorization.
int verbose_
Toggles the output level.
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
void PrintLine() const
Prints line on std::cout.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().