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 ) ) ;
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");
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.