31#include "Epetra_MpiComm.h"
33 This code cannot be compiled without mpi.h.
34#include
"Epetra_Comm.h"
36#include "Epetra_Map.h"
37#include "Epetra_LocalMap.h"
38#include "Epetra_Vector.h"
39#include "Epetra_RowMatrix.h"
40#include "Epetra_Operator.h"
41#include "Epetra_Import.h"
42#include "Epetra_Export.h"
43#include "Epetra_CrsMatrix.h"
44#include "supermatrix.h"
45#include "superlu_ddefs.h"
50#include "Comm_assert_equal.h"
53#ifdef EPETRA_CRSMATRIX_CONSTRUCT_FROM_ROWMATRIX
54#include "ExtractCrsFromRowMatrix.h"
69 for ( i = 1; i*i <= NumProcs; i++ )
72 for ( numrows = i-1 ; done == false ; ) {
73 int NumCols = NumProcs / numrows ;
74 if ( numrows * NumCols == NumProcs )
109 SUPERLU_FREE(
A.Store);
176 bool CheckExtraction =
false;
178 Epetra_RowMatrix *RowMatrixA =
179 dynamic_cast<Epetra_RowMatrix *
>(
Problem_->GetOperator());
183 Epetra_CrsMatrix *CastCrsMatrixA =
dynamic_cast<Epetra_CrsMatrix*
>(RowMatrixA) ;
184#ifdef EPETRA_CRSMATRIX_CONSTRUCT_FROM_ROWMATRIX
185 Epetra_CrsMatrix *ExtractCrsMatrixA = 0;
187 Epetra_CrsMatrix *Phase2Mat = 0 ;
188 const Epetra_Comm &Comm = RowMatrixA->Comm();
190 int iam = Comm.MyPID() ;
196 if ( iam == 0 ) cin >> hatever ;
206 if ( CastCrsMatrixA != 0 && ! CheckExtraction ) {
207 Phase2Mat = CastCrsMatrixA ;
209#ifdef EPETRA_CRSMATRIX_CONSTRUCT_FROM_ROWMATRIX
210 ExtractCrsMatrixA =
new Epetra_CrsMatrix( *RowMatrixA ) ;
212 Phase2Mat = ExtractCrsMatrixA ;
214 if ( CheckExtraction )
215 assert( CrsMatricesAreIdentical( CastCrsMatrixA, ExtractCrsMatrixA ) ) ;
222 assert( Phase2Mat != NULL ) ;
223 const Epetra_Map &Phase2Matmap = Phase2Mat->RowMap() ;
228 const Epetra_MpiComm & comm1 =
dynamic_cast<const Epetra_MpiComm &
> (Comm);
229 MPI_Comm MPIC = comm1.Comm() ;
231 int IsLocal = ( Phase2Matmap.NumMyElements() ==
232 Phase2Matmap.NumGlobalElements() )?1:0;
233 Comm.Broadcast( &IsLocal, 1, 0 ) ;
235 assert( Comm_assert_equal( &Comm, IsLocal ) );
238 Epetra_CrsMatrix *Phase3Mat = 0 ;
240 int NumGlobalElements_ = Phase2Matmap.NumGlobalElements() ;
244 int NumMyElements_ = 0 ;
245 if (iam==0) NumMyElements_ = NumGlobalElements_;
246 Epetra_Map SerialMap( NumGlobalElements_, NumMyElements_, 0, Comm );
247 Epetra_CrsMatrix SerialCrsMatrixA(Copy, SerialMap, 0);
251 Phase3Mat = Phase2Mat ;
254 Epetra_Export export_to_serial( Phase2Matmap, SerialMap);
256 SerialCrsMatrixA.Export( *Phase2Mat, export_to_serial, Add );
258 SerialCrsMatrixA.FillComplete() ;
259 Phase3Mat = &SerialCrsMatrixA ;
269 const Epetra_Map &Phase3Matmap = Phase3Mat->RowMap() ;
271 numrows = Phase3Mat->NumGlobalRows();
272 int numcols = Phase3Mat->NumGlobalCols();
273 int numentries = Phase3Mat->NumGlobalNonzeros();
275 Epetra_CrsMatrix Phase3MatTrans(Copy, Phase3Matmap, 0);
276 Epetra_CrsMatrix *Phase4Mat;
279 Phase4Mat = Phase3Mat ;
282 Phase4Mat = &Phase3MatTrans ;
289 int * AllIDs =
new int[
numrows];
290 for (
int i = 0; i <
numrows ; i++ ) AllIDs[i] = i ;
293 Epetra_Map ReplicatedMap( -1,
numrows, AllIDs, 0, Comm);
298 Epetra_Import importer( ReplicatedMap, Phase3Matmap );
300 Epetra_CrsMatrix Phase5Mat(Copy, ReplicatedMap, 0);
301 EPETRA_CHK_ERR( Phase5Mat.Import( *Phase4Mat, importer, Insert) );
304 assert( Phase5Mat.NumMyRows() == Phase4Mat->NumGlobalRows() ) ;
309 Epetra_MultiVector *vecX =
Problem_->GetLHS() ;
310 Epetra_MultiVector *vecB =
Problem_->GetRHS() ;
317 nrhs = vecX->NumVectors() ;
321 int nArows = Phase3Mat->NumGlobalRows() ;
322 int nAcols = Phase3Mat->NumGlobalCols() ;
324#ifdef ONE_VECTOR_ONLY
325 assert( vecX->NumVectors() == 1 ) ;
326 assert( vecB->NumVectors() == 1 ) ;
328 Epetra_Vector *vecXvector =
dynamic_cast<Epetra_Vector*
>(vecX) ;
329 Epetra_Vector *vecBvector =
dynamic_cast<Epetra_Vector*
>(vecB) ;
331 assert( vecXvector != 0 ) ;
332 assert( vecBvector != 0 ) ;
334 Epetra_Vector vecXreplicated( ReplicatedMap ) ;
335 Epetra_Vector vecBreplicated( ReplicatedMap ) ;
337 Epetra_MultiVector *vecXvector = (vecX) ;
338 Epetra_MultiVector *vecBvector = (vecB) ;
340 Epetra_MultiVector vecXreplicated( ReplicatedMap, nrhs ) ;
341 Epetra_MultiVector vecBreplicated( ReplicatedMap, nrhs ) ;
344 Epetra_Import ImportToReplicated( ReplicatedMap, Phase2Matmap);
346 vecXreplicated.Import( *vecXvector, ImportToReplicated, Insert ) ;
347 vecBreplicated.Import( *vecBvector, ImportToReplicated, Insert ) ;
349 assert( nArows == vecXreplicated.MyLength() ) ;
350 assert( nAcols == vecBreplicated.MyLength() ) ;
356 assert( vecBreplicated.ExtractView( &bValues, &bLda ) == 0 ) ;
357 assert( vecXreplicated.ExtractView( &xValues, &xLda ) == 0 ) ;
364 Ai.resize( EPETRA_MAX(
numrows, numentries) ) ;
372 for ( MyRow = 0; MyRow <
numrows; MyRow++ ) {
373 int status = Phase5Mat.ExtractMyRowView( MyRow, NumEntries, RowValues, ColIndices ) ;
374 assert( status == 0 ) ;
375 Ap[MyRow] = Ai_index ;
376 for (
int j = 0; j < NumEntries; j++ ) {
377 Ai[Ai_index] = ColIndices[j] ;
378 Aval[Ai_index] = RowValues[j] ;
401 assert( Comm_assert_equal( &Comm, numentries ) );
402 assert( Comm_assert_equal( &Comm,
numrows ) );
403 assert( Comm_assert_equal( &Comm, numcols ) );
407 assert(
numprocs == Comm.NumProc() ) ;
419 for (
int ii = 0; ii < min( numentries, 10 ) ; ii++ ) {
420 assert( Comm_assert_equal( &Comm,
Aval[ii] ) ) ;
422 for (
int ii = 0; ii < min( numcols+1, 10 ) ; ii++ ) {
423 assert( Comm_assert_equal( &Comm,
Ai[ii] ) );
424 assert( Comm_assert_equal( &Comm,
Ap[ii] ) );
426 for (
int ii = 0; ii < min(
numrows, 10 ) ; ii++ ) {
427 assert( Comm_assert_equal( &Comm, bValues[ii] ) );
434 if ( !(
berr = doubleMalloc_dist(nrhs)) )
438 dCreate_CompCol_Matrix_dist(&
A,
numrows, numcols,
439 numentries, &
Aval[0], &
Ai[0],
440 &
Ap[0], SLU_NC, SLU_D, SLU_GE);
444 cout <<
" Here is A " << endl ;
445 dPrint_CompCol_Matrix_dist( &
A );
446 cout <<
" That was A " <<
"numrows = " <<
numrows << endl ;
447 cout <<
"numcols = " << numcols << endl ;
454 assert(
options.Fact == DOFACT );
466 for (
int j = 0 ; j < nrhs; j++ )
467 for (
int i = 0 ; i <
numrows; i++ ) xValues[i+j*xLda] = bValues[i+j*xLda];
488 for (
int i = 0 ; i <
numrows; i++ ) {
490 lid[0] = Phase2Matmap.LID( i ) ;
492 for (
int j = 0 ; j < nrhs; j++ ) {
493 vecXvector->ReplaceMyValue( lid[0], j, xValues[i + ldb * j] ) ;
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
#define EPETRA_CHK_ERR(xxx)
int SLU_NumRows(int NumProcs)
SuperludistOO(const Epetra_LinearProblem &LinearProblem)
SuperludistOO Constructor.
virtual ~SuperludistOO(void)
SuperludistOO Destructor.
int Solve(bool Factor)
All computation is performed during the call to Solve()
bool GetTrans() const
Return the transpose flag.
ScalePermstruct_t ScalePermstruct
const Epetra_LinearProblem * Problem_
superlu_options_t options