31 #include "Epetra_MpiComm.h" 33 #include "Epetra_Comm.h" 35 #include "Epetra_Map.h" 36 #include "Epetra_Vector.h" 37 #include "Epetra_RowMatrix.h" 38 #include "Epetra_Operator.h" 39 #include "Epetra_Import.h" 40 #include "Epetra_CrsMatrix.h" 48 #include "BridgeMPI.h" 53 Epetra_MultiVector * X,
54 Epetra_MultiVector *
B) {
129 FILE *msgFile = fopen(
"msgFile",
"w" ) ;
133 const Epetra_Comm &Comm = RowMatrixA->Comm();
134 int iam = Comm.MyPID() ;
137 if (
verbose ) matFile = fopen(
"SPO.m",
"w" ) ;
139 Epetra_CrsMatrix *matA =
dynamic_cast<Epetra_CrsMatrix*
>(RowMatrixA) ;
140 assert( matA !=
NULL ) ;
143 Epetra_MultiVector *vecX =
GetLHS() ;
144 Epetra_MultiVector *vecB =
GetRHS() ;
146 const Epetra_Map &matAmap = matA->RowMap() ;
147 const Epetra_MpiComm & comm1 =
dynamic_cast<const Epetra_MpiComm &
> (Comm);
148 MPI_Comm MPIC = comm1.Comm() ;
150 int IsLocal = ( matAmap.NumMyElements() ==
151 matAmap.NumGlobalElements() )?1:0;
152 Comm.Broadcast( &IsLocal, 1, 0 ) ;
157 int nArows = matA->NumGlobalRows() ;
158 int nAcols = matA->NumGlobalCols() ;
159 int nproc = Comm.NumProc() ;
161 assert( vecX->NumVectors() == 1 ) ;
162 assert( vecB->NumVectors() == 1 ) ;
169 assert( matAmap.NumMyElements() == matAmap.NumGlobalElements() ) ;
171 assert( matAmap.NumMyElements() == 0 ) ;
174 int numentries = matA->NumGlobalNonzeros();
175 vector <int>rowindices( numentries ) ;
176 vector <int>colindices( numentries ) ;
177 vector <double>values( numentries ) ;
182 int numrows = matA->NumGlobalRows();
198 for (
int MyRow = 0; MyRow <numrows; MyRow++ ) {
199 assert( matA->ExtractMyRowView( MyRow, NumEntries, RowValues, ColIndices ) == 0 ) ;
200 for (
int j = 0; j < NumEntries; j++ ) {
201 rowindices[entrynum] = MyRow ;
202 colindices[entrynum] = ColIndices[j] ;
203 values[entrynum] = RowValues[j] ;
207 mtxA = InpMtx_new() ;
209 InpMtx_init( mtxA, 1, SPOOLES_REAL, 0, 0 ) ;
213 InpMtx_inputRealTriples( mtxA, numentries, &colindices[0],
214 &rowindices[0], &values[0] ) ;
216 InpMtx_inputRealTriples( mtxA, numentries, &rowindices[0],
217 &colindices[0], &values[0] ) ;
235 assert( vecB->ExtractView( &bValues, &bLda ) == 0 ) ;
236 assert( vecX->ExtractView( &xValues, &xLda ) == 0 ) ;
242 mtxX = DenseMtx_new() ;
243 mtxY = DenseMtx_new() ;
245 DenseMtx_init( mtxX, SPOOLES_REAL, -1, -1, numrows, 1, 1, numrows );
246 DenseMtx_init( mtxY, SPOOLES_REAL, -1, -1, numrows, 1, 1, numrows );
248 DenseMtx_init( mtxX, SPOOLES_REAL, 0, 0, numrows, 1, 1, numrows );
249 DenseMtx_init( mtxY, SPOOLES_REAL, 0, 0, numrows, 1, 1, numrows );
252 DenseMtx_dimensions( mtxY, &nrow, &nnrhs ) ;
253 assert( nrow = numrows ) ;
254 assert( nnrhs == 1 ) ;
257 if (
verbose) InpMtx_writeForMatlab( mtxA,
"mtxA", matFile ) ;
262 for (
int i = 0 ; i < numrows; i ++ )
264 DenseMtx_setRealEntry( mtxY, i, 0, bValues[i] );
266 if (
verbose ) DenseMtx_writeForMatlab( mtxY,
"mtxY", matFile ) ;
267 DenseMtx_zero(mtxX) ;
271 int type = SPOOLES_REAL ;
272 int symmetryflag = SPOOLES_NONSYMMETRIC ;
282 BridgeMPI *bridgeMPI = BridgeMPI_new() ;
283 BridgeMPI_setMPIparams(bridgeMPI, nproc, iam, MPIC ) ;
284 BridgeMPI_setMatrixParams(bridgeMPI, numrows, type, symmetryflag) ;
286 double droptol = 0.0 ;
288 PatchAndGoInfo *PatchAndGo = 0 ;
289 BridgeMPI_setFactorParams(bridgeMPI, FRONTMTX_DENSE_FRONTS,
290 SPOOLES_PIVOTING, tau, droptol, lookahead,
292 BridgeMPI_setMessageInfo(bridgeMPI, msglvl, msgFile) ;
293 assert( type == SPOOLES_REAL ) ;
294 assert( msglvl >= 0 && msglvl <= 2 ) ;
296 rc = BridgeMPI_setup(bridgeMPI, mtxA) ;
298 fprintf(stderr,
"\n error return %d from BridgeMPI_setup()", rc) ;
303 int nfront, nfind, nfent, nsolveops;
305 rc = BridgeMPI_factorStats(bridgeMPI, type, symmetryflag, &nfront,
306 &nfind, &nfent, &nsolveops, &nfactorops) ;
309 "\n error return %d from BridgeMPI_factorStats()", rc) ;
318 rc = BridgeMPI_factorSetup(bridgeMPI, 0, 0.0) ;
320 std::stringstream Message ;
322 Message <<
" SPOOLES factorsetup failed with return code " <<
324 string mess = Message.str() ;
333 int permuteflag = 1 ;
335 rc = BridgeMPI_factor(bridgeMPI, mtxA, permuteflag, &errorcode) ;
336 assert( permuteflag == 1 ) ;
339 std::stringstream Message ;
341 Message <<
" SPOOLES factorization failed with return code " <<
342 rc <<
" and error code " << errorcode ;
343 string mess = Message.str() ;
352 rc = BridgeMPI_solveSetup(bridgeMPI) ;
354 std::stringstream Message ;
356 Message <<
" SPOOLES BridgeMPI_solveSetup failed with return code" <<
358 string mess = Message.str() ;
362 assert( permuteflag == 1 ) ;
363 rc = BridgeMPI_solve(bridgeMPI, permuteflag, mtxX, mtxY) ;
364 assert( permuteflag == 1 ) ;
366 std::stringstream Message ;
368 Message <<
" SPOOLES BridgeMPI_solve failed with return code" <<
370 string mess = Message.str() ;
374 if (
verbose ) DenseMtx_writeForMatlab( mtxX,
"mtxXX", matFile ) ;
375 if (
verbose ) fclose(matFile);
382 for (
int i = 0 ; i < numrows; i ++ )
384 DenseMtx_realEntry( mtxX, i, 0, &xValues[i] );
392 DenseMtx_free(mtxX) ;
393 DenseMtx_free(mtxY) ;
395 BridgeMPI_free(bridgeMPI) ;
int SetRHS(Epetra_MultiVector *B)
Epetra_RowMatrix * UserMatrix_
Epetra_RowMatrix * GetUserMatrix() const
int SetUserMatrix(Epetra_RowMatrix *UserMatrix)
#define EPETRA_CHK_ERR(xxx)
int SetLHS(Epetra_MultiVector *X)
Epetra_Operator * UserOperator_
Epetra_MultiVector * GetLHS() const
Epetra_MultiVector * GetRHS() const