32 #include "Epetra_MpiComm.h" 34 This code cannot be compiled without mpi.h.
35 #include
"Epetra_Comm.h" 37 #include "Epetra_Map.h" 38 #include "Epetra_LocalMap.h" 39 #include "Epetra_Vector.h" 40 #include "Epetra_RowMatrix.h" 41 #include "Epetra_Operator.h" 42 #include "Epetra_Import.h" 43 #include "Epetra_Export.h" 44 #include "Epetra_CrsMatrix.h" 49 #include "Comm_assert_equal.h" 52 #ifdef EPETRA_CRSMATRIX_CONSTRUCT_FROM_ROWMATRIX 53 #include "ExtractCrsFromRowMatrix.h" 68 for ( i = 1; i*i <= NumProcs; i++ )
71 for ( NumProcRows = i-1 ; done == false ; ) {
72 int NumCols = NumProcs / NumProcRows ;
73 if ( NumProcRows * NumCols == NumProcs )
108 SUPERLU_FREE(
A.Store );
112 if (
options.SolveInitialized ) {
115 superlu_gridexit(&
grid);
173 #undef EPETRA_CHK_ERR 174 #define EPETRA_CHK_ERR(xxx) assert( (xxx) == 0 ) 183 bool CheckExtraction =
false;
187 Epetra_RowMatrix *RowMatrixA =
188 dynamic_cast<Epetra_RowMatrix *
>(
Problem_->GetOperator());
192 Epetra_CrsMatrix *CastCrsMatrixA =
dynamic_cast<Epetra_CrsMatrix*
>(RowMatrixA) ;
193 #ifdef EPETRA_CRSMATRIX_CONSTRUCT_FROM_ROWMATRIX 194 Epetra_CrsMatrix *ExtractCrsMatrixA = 0;
196 Epetra_CrsMatrix *Phase2Mat = 0 ;
197 const Epetra_Comm &Comm = RowMatrixA->Comm();
199 int iam = Comm.MyPID() ;
205 if ( iam == 0 ) cin >> hatever ;
216 if ( CastCrsMatrixA != 0 && ! CheckExtraction ) {
217 Phase2Mat = CastCrsMatrixA ;
219 #ifdef EPETRA_CRSMATRIX_CONSTRUCT_FROM_ROWMATRIX 220 ExtractCrsMatrixA =
new Epetra_CrsMatrix( *RowMatrixA ) ;
222 Phase2Mat = ExtractCrsMatrixA ;
224 if ( CheckExtraction )
225 assert( CrsMatricesAreIdentical( CastCrsMatrixA, ExtractCrsMatrixA ) ) ;
231 assert( Phase2Mat !=
NULL ) ;
232 const Epetra_Map &Phase2Matmap = Phase2Mat->RowMap() ;
247 cout <<
" A Comm.NumProc() = " << Comm.NumProc() << endl ;
249 cout <<
" true = " <<
true <<
" LInMap = " << Phase2Matmap.LinearMap() << endl ;
250 cout <<
" Superludist2_OO.cpp:: traceback mode = " << Epetra_Object::GetTracebackMode() << endl ;
251 cerr <<
" Send this to cerr cerr cerr traceback mode = " << Epetra_Object::GetTracebackMode() << endl ;
254 numrows = Phase2Matmap.NumGlobalElements() ;
255 assert(
numrows == Phase2Mat->NumGlobalRows() ) ;
256 int numcols = Phase2Mat->NumGlobalCols() ;
264 int m_per_p =
numrows / Comm.NumProc() ;
265 int remainder =
numrows - ( m_per_p * Comm.NumProc() );
266 int MyFirstElement = iam * m_per_p + EPETRA_MIN( iam, remainder );
267 int MyFirstNonElement = (iam+1) * m_per_p + EPETRA_MIN( iam+1, remainder );
268 int NumExpectedElements = MyFirstNonElement - MyFirstElement ;
271 int IsLocal = ( Phase2Matmap.NumMyElements() ==
272 Phase2Matmap.NumGlobalElements() )?1:0;
273 Comm.Broadcast( &IsLocal, 1, 0 ) ;
277 Epetra_CrsMatrix *Phase3Mat = 0 ;
278 Epetra_Map DistMap( Phase2Matmap.NumGlobalElements(), NumExpectedElements, 0, Comm );
279 Epetra_CrsMatrix DistCrsMatrixA(
Copy, DistMap, 0);
281 bool redistribute = true ;
282 if ( redistribute ) {
284 Epetra_Export export_to_dist( Phase2Matmap, DistMap);
286 DistCrsMatrixA.Export( *Phase2Mat, export_to_dist, Add );
288 DistCrsMatrixA.FillComplete() ;
289 Phase3Mat = &DistCrsMatrixA ;
299 int Phase2NumElements = Phase2Matmap.NumMyElements() ;
300 vector <int> MyRows( Phase2NumElements ) ;
301 Phase2Matmap.MyGlobalElements( &MyRows[0] ) ;
302 for (
int row = 0 ; row < Phase2NumElements ; row++ ) {
306 Phase3Mat = Phase2Mat ;
310 assert( MyFirstElement == Phase3Mat->RowMap().MinMyGID() ) ;
311 assert( NumExpectedElements == Phase3Mat->RowMap().NumMyElements() ) ;
312 assert( MyFirstElement+NumExpectedElements-1 == Phase3Mat->RowMap().MaxMyGID() ) ;
315 int MyActualFirstElement = Phase3Mat->RowMap().MinMyGID() ;
316 int NumMyElements = Phase3Mat->NumMyRows() ;
317 vector <int> MyRowPtr( NumMyElements+1 ) ;
322 int CurrentRowPtr = 0 ;
323 for (
int i = 0; i < NumMyElements ; i++ ) {
324 CurrentRowPtr += Phase3Mat->NumMyEntries( i ) ;
325 MyRowPtr[i+1] = CurrentRowPtr ;
332 int nnz_loc = Phase3Mat->NumMyNonzeros() ;
337 Ap.resize( NumMyElements+1 );
338 Ai.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ;
339 Aval.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ;
346 int num_my_cols = Phase3Mat->NumMyCols() ;
347 vector <int>Global_Columns( num_my_cols ) ;
348 for (
int i = 0 ; i < num_my_cols ; i ++ ) {
349 Global_Columns[i] = Phase3Mat->GCID( i ) ;
352 for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
353 int status = Phase3Mat->ExtractMyRowView( MyRow, NzThisRow, RowValues, ColIndices ) ;
354 assert( status == 0 ) ;
355 Ap[MyRow] = Ai_index ;
356 assert(
Ap[MyRow] == MyRowPtr[MyRow] ) ;
357 for (
int j = 0; j < NzThisRow; j++ ) {
358 Ai[Ai_index] = Global_Columns[ColIndices[j]] ;
359 Aval[Ai_index] = RowValues[j] ;
363 assert( NumMyElements == MyRow );
364 Ap[ NumMyElements ] = Ai_index ;
371 Epetra_MultiVector *vecX =
Problem_->GetLHS() ;
372 Epetra_MultiVector *vecB =
Problem_->GetRHS() ;
379 nrhs = vecX->NumVectors() ;
383 Epetra_MultiVector vecXdistributed( DistMap, nrhs ) ;
384 Epetra_MultiVector vecBdistributed( DistMap, nrhs ) ;
391 Epetra_MultiVector* vecXptr;
392 Epetra_MultiVector* vecBptr;
394 if ( redistribute ) {
395 Epetra_Import ImportToDistributed( DistMap, Phase2Matmap);
397 vecXdistributed.Import( *vecX, ImportToDistributed, Insert ) ;
398 vecBdistributed.Import( *vecB, ImportToDistributed, Insert ) ;
400 vecXptr = &vecXdistributed ;
401 vecBptr = &vecBdistributed ;
420 const Epetra_MpiComm & comm1 =
dynamic_cast<const Epetra_MpiComm &
> (Comm);
421 MPI_Comm MPIC = comm1.Comm() ;
432 assert(
numprocs == Comm.NumProc() ) ;
441 dCreate_CompRowLoc_Matrix_dist( &
A,
numrows, numcols,
442 nnz_loc, NumMyElements, MyActualFirstElement,
444 SLU_NR_loc, SLU_D, SLU_GE );
452 assert(
options.Fact == DOFACT );
464 for (
int j = 0 ; j < nrhs; j++ )
465 for (
int i = 0 ; i < NumMyElements; i++ ) xValues[i+j*ldx] = bValues[i+j*ldb];
470 vector<double>berr(nrhs);
479 if ( redistribute ) {
480 Epetra_Import ImportBackToOriginal( Phase2Matmap,DistMap);
482 vecX->Import( *vecXptr, ImportBackToOriginal, Insert ) ;
SOLVEstruct_t SOLVEstruct
virtual ~Superludist2_OO(void)
Superludist2_OO Destructor.
#define EPETRA_CHK_ERR(xxx)
superlu_options_t options
bool GetTrans() const
Return the transpose flag.
ScalePermstruct_t ScalePermstruct
Superludist2_OO(const Epetra_LinearProblem &LinearProblem)
Superludist2_OO Constructor.
const Epetra_LinearProblem * Problem_
int Solve(bool Factor)
All computation is performed during the call to Solve()
int SLU_NumProcRows(int NumProcs)