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 "Teuchos_RCP.hpp"
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();
125 Teuchos::RCP<Epetra_CrsGraph> Graph;
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;
329 Epetra_CrsMatrix DscMat(Copy,
DscRowMap(), 0);
333 DscColMap_ = Teuchos::rcp(
new Epetra_Map(DscMat.RowMatrixColMap()));
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();
558 int nnz =
GetProblem()->GetMatrix()->NumGlobalNonzeros();
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;
#define AMESOS_CHK_ERR(a)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int PerformSymbolicFactorization()
Performs the symbolic factorization.
int PerformNumericFactorization()
Performs the numeric factorization.
bool MatrixShapeOK() const
Returns true if DSCPACK can handle this matrix shape.
RCP< Epetra_Map > DscRowMap_
bool A_and_LU_built
Tells us whether to free them.
void PrintTiming() const
Prints timing information.
RCP< Epetra_Import > Importer_
const Epetra_Import & Importer() const
const Epetra_Map & DscColMap() const
Teuchos::RCP< Amesos_Dscpack_Pimpl > PrivateDscpackData_
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
int Solve()
Solves A X = B (or AT x = B)
void PrintStatus() const
Prints information about the factorization and solution phases.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
RCP< Epetra_Map > DscColMap_
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Amesos_Dscpack(const Epetra_LinearProblem &LinearProblem)
Amesos_Dscpack Constructor.
const Epetra_Map & DscRowMap() const
~Amesos_Dscpack(void)
Amesos_Dscpack Destructor.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
const Epetra_LinearProblem * Problem_
Pointer to the linear problem.
int * GlobalStructNewColNum
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.
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.