Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
SpoolesOO.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: Direct Sparse Solver Package
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25//
26// ***********************************************************************
27// @HEADER
28
29#include "SpoolesOO.h"
30#ifdef EPETRA_MPI
31#include "Epetra_MpiComm.h"
32#else
33#include "Epetra_Comm.h"
34#endif
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"
41#include <sstream>
42#include <vector>
43
44//
45// SPOOLES include file:
46//
47extern "C" {
48#include "BridgeMPI.h"
49}
50
51//=============================================================================
52SpoolesOO::SpoolesOO(Epetra_RowMatrix * A,
53 Epetra_MultiVector * X,
54 Epetra_MultiVector * B) {
55 // AllocAzArrays();
57
58 inConstructor_ = true; // Shut down complaints about zero pointers for a while
60
61 SetLHS(X);
62 SetRHS(B);
63 inConstructor_ = false;
64}
65
66//=============================================================================
68 // AllocAzArrays();
70}
71
72//=============================================================================
74 // DeleteMemory();
75 // DeleteAzArrays();
76}
77
78//=============================================================================
79int SpoolesOO::SetUserMatrix(Epetra_RowMatrix * UserMatrix) {
80
81 if (UserMatrix == 0 && inConstructor_ == true) return(0);
82 if (UserMatrix == 0) EPETRA_CHK_ERR(-1);
83
84 UserMatrix_ = UserMatrix;
85
86 return(0);
87}
88
89//=============================================================================
90int SpoolesOO::SetLHS(Epetra_MultiVector * X) {
91
92 if (X == 0 && inConstructor_ == true) return(0);
93 if (X == 0) EPETRA_CHK_ERR(-1);
94 X_ = X;
95 X_->ExtractView(&x_, &x_LDA_);
96 return(0);
97}
98//=============================================================================
99int SpoolesOO::SetRHS(Epetra_MultiVector * B) {
100
101 if (B == 0 && inConstructor_ == true) return(0);
102 if (B == 0) EPETRA_CHK_ERR(-1);
103 B_ = B;
104 B_->ExtractView(&b_, &b_LDA_);
105
106 return(0);
107}
109
110 UserOperator_ = 0;
111 UserMatrix_ = 0;
112 // PrecOperator_ = 0;
113 // PrecMatrix_ = 0;
114 X_ = 0;
115 B_ = 0;
116
117 x_LDA_ = 0;
118 x_ = 0;
119 b_LDA_ = 0;
120 b_ = 0;
121
122 return(0);
123
124}
125
126//=============================================================================
127
129 FILE *msgFile = fopen("msgFile", "w" ) ;
130 FILE *matFile = 0 ;
131
132 Epetra_RowMatrix *RowMatrixA = (GetUserMatrix()) ;
133 const Epetra_Comm &Comm = RowMatrixA->Comm();
134 int iam = Comm.MyPID() ;
135 bool verbose = false ; // Other option is (iam == 0 ) ;
136 verbose = ( iam == 0 ) ;
137 if ( verbose ) matFile = fopen("SPO.m", "w" ) ;
138
139 Epetra_CrsMatrix *matA = dynamic_cast<Epetra_CrsMatrix*>(RowMatrixA) ;
140 assert( matA != NULL ) ; // SuperludistOO shows how to convert a
141 // non-CrsMatrix into a CrsMatrix
142
143 Epetra_MultiVector *vecX = GetLHS() ;
144 Epetra_MultiVector *vecB = GetRHS() ;
145
146 const Epetra_Map &matAmap = matA->RowMap() ;
147 const Epetra_MpiComm & comm1 = dynamic_cast<const Epetra_MpiComm &> (Comm);
148 MPI_Comm MPIC = comm1.Comm() ;
149
150 int IsLocal = ( matAmap.NumMyElements() ==
151 matAmap.NumGlobalElements() )?1:0;
152 Comm.Broadcast( &IsLocal, 1, 0 ) ;
153
154 EPETRA_CHK_ERR( 1 - IsLocal ); // SuperludistOO shows hows to
155 // deal with distributed matrices.
156
157 int nArows = matA->NumGlobalRows() ;
158 int nAcols = matA->NumGlobalCols() ;
159 int nproc = Comm.NumProc() ;
160
161 assert( vecX->NumVectors() == 1 ) ;
162 assert( vecB->NumVectors() == 1 ) ;
163
164 // assert( nArows == vecX->MyLength() ) ;
165 // assert( nAcols == vecB->MyLength() ) ;
166
167 // assert( matAmap.DistributedGlobal() == false ) ;
168 if ( iam == 0 ) {
169 assert( matAmap.NumMyElements() == matAmap.NumGlobalElements() ) ;
170 } else {
171 assert( matAmap.NumMyElements() == 0 ) ;
172 }
173
174 int numentries = matA->NumGlobalNonzeros();
175 vector <int>rowindices( numentries ) ;
176 vector <int>colindices( numentries ) ;
177 vector <double>values( numentries ) ;
178 int NumEntries;
179 double *RowValues;
180
181 int *ColIndices;
182 int numrows = matA->NumGlobalRows();
183 int entrynum = 0 ;
184
185 //
186 // The following lines allow us time to attach the debugger
187 //
188 char hatever;
189 // if ( iam == 0 ) cin >> hatever ;
190 Comm.Barrier();
191
192
193 InpMtx *mtxA;
194 if ( iam==0 ) {
195 //
196 // Create mtxA on proc 0
197 //
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] ;
204 entrynum++;
205 }
206 }
207 mtxA = InpMtx_new() ;
208
209 InpMtx_init( mtxA, 1, SPOOLES_REAL, 0, 0 ) ; // I can't figure out
210 // what the right bound for the number of vectors is - Ken
211
212 if ( GetTrans() ) {
213 InpMtx_inputRealTriples( mtxA, numentries, &colindices[0],
214 &rowindices[0], &values[0] ) ;
215 } else {
216 InpMtx_inputRealTriples( mtxA, numentries, &rowindices[0],
217 &colindices[0], &values[0] ) ;
218 }
219 } else {
220 mtxA = 0 ;
221 }
222
223 // return 1; OK to here
224
225 DenseMtx *mtxX = 0 ;
226 DenseMtx *mtxY = 0 ;
227 double *bValues ;
228 double *xValues ;
229 if ( iam == 0 ) {
230 //
231 // Convert Epetra_Vector x and Epetra_Vector b arrays
232 //
233 int bLda, xLda ;
234
235 assert( vecB->ExtractView( &bValues, &bLda ) == 0 ) ;
236 assert( vecX->ExtractView( &xValues, &xLda ) == 0 ) ;
237
238
239 //
240 // SPOOLES matrices
241 //
242 mtxX = DenseMtx_new() ;
243 mtxY = DenseMtx_new() ;
244#ifdef OLDWAY
245 DenseMtx_init( mtxX, SPOOLES_REAL, -1, -1, numrows, 1, 1, numrows );
246 DenseMtx_init( mtxY, SPOOLES_REAL, -1, -1, numrows, 1, 1, numrows );
247#else
248 DenseMtx_init( mtxX, SPOOLES_REAL, 0, 0, numrows, 1, 1, numrows );
249 DenseMtx_init( mtxY, SPOOLES_REAL, 0, 0, numrows, 1, 1, numrows );
250#endif
251 int nrow, nnrhs ;
252 DenseMtx_dimensions( mtxY, &nrow, &nnrhs ) ;
253 assert( nrow = numrows ) ;
254 assert( nnrhs == 1 ) ;
255
256
257 if (verbose) InpMtx_writeForMatlab( mtxA, "mtxA", matFile ) ;
258 //
259 // This is a maximally inefficient way to create the right hand side
260 // matrix, but I could not find anything better in the documentation.
261 //
262 for ( int i = 0 ; i < numrows; i ++ )
263 {
264 DenseMtx_setRealEntry( mtxY, i, 0, bValues[i] );
265 }
266 if ( verbose ) DenseMtx_writeForMatlab( mtxY, "mtxY", matFile ) ;
267 DenseMtx_zero(mtxX) ;
268 }
269
270
271 int type = SPOOLES_REAL ;
272 int symmetryflag = SPOOLES_NONSYMMETRIC ;
273 // SPOOLES requires a message level and a message File
274 int msglvl = 0;
275 int rc; // return code
276
277 /*
278 --------------------------------
279 create and setup a Bridge object
280 --------------------------------
281 */
282 BridgeMPI *bridgeMPI = BridgeMPI_new() ;
283 BridgeMPI_setMPIparams(bridgeMPI, nproc, iam, MPIC ) ;
284 BridgeMPI_setMatrixParams(bridgeMPI, numrows, type, symmetryflag) ;
285 double tau = 100.0 ;
286 double droptol = 0.0 ;
287 int lookahead = 0 ;
288 PatchAndGoInfo *PatchAndGo = 0 ;
289 BridgeMPI_setFactorParams(bridgeMPI, FRONTMTX_DENSE_FRONTS,
290 SPOOLES_PIVOTING, tau, droptol, lookahead,
291 PatchAndGo ) ;
292 BridgeMPI_setMessageInfo(bridgeMPI, msglvl, msgFile) ;
293 assert( type == SPOOLES_REAL ) ;
294 assert( msglvl >= 0 && msglvl <= 2 ) ;
295
296 rc = BridgeMPI_setup(bridgeMPI, mtxA) ;
297 if ( rc != 1 ) {
298 fprintf(stderr, "\n error return %d from BridgeMPI_setup()", rc) ;
299 exit(-1) ;
300 }
301 Comm.Barrier();
302
303 int nfront, nfind, nfent, nsolveops;
304 double nfactorops;
305 rc = BridgeMPI_factorStats(bridgeMPI, type, symmetryflag, &nfront,
306 &nfind, &nfent, &nsolveops, &nfactorops) ;
307 if ( rc != 1 ) {
308 fprintf(stderr,
309 "\n error return %d from BridgeMPI_factorStats()", rc) ;
310 exit(-1) ;
311 }
312
313 /*
314 --------------------------------
315 setup the parallel factorization
316 --------------------------------
317 */
318 rc = BridgeMPI_factorSetup(bridgeMPI, 0, 0.0) ;
319 if (rc != 1 ) {
320 std::stringstream Message ;
321
322 Message << " SPOOLES factorsetup failed with return code " <<
323 rc ;
324 string mess = Message.str() ;
325 throw( mess ) ;
326 }
327
328 /*
329 -----------------
330 factor the matrix
331 -----------------
332 */
333 int permuteflag = 1 ;
334 int errorcode ;
335 rc = BridgeMPI_factor(bridgeMPI, mtxA, permuteflag, &errorcode) ;
336 assert( permuteflag == 1 ) ;
337
338 if (rc != 1 ) {
339 std::stringstream Message ;
340
341 Message << " SPOOLES factorization failed with return code " <<
342 rc << " and error code " << errorcode ;
343 string mess = Message.str() ;
344 throw( mess ) ;
345 }
346 /*
347 ----------------
348 solve the system
349 ----------------
350 */
351
352 rc = BridgeMPI_solveSetup(bridgeMPI) ;
353 if (rc != 1 ) {
354 std::stringstream Message ;
355
356 Message << " SPOOLES BridgeMPI_solveSetup failed with return code" <<
357 rc << endl ;
358 string mess = Message.str() ;
359 throw( mess ) ;
360 }
361
362 assert( permuteflag == 1 ) ;
363 rc = BridgeMPI_solve(bridgeMPI, permuteflag, mtxX, mtxY) ;
364 assert( permuteflag == 1 ) ;
365 if (rc != 1 ) {
366 std::stringstream Message ;
367
368 Message << " SPOOLES BridgeMPI_solve failed with return code" <<
369 rc << endl ;
370 string mess = Message.str() ;
371 throw( mess ) ;
372 }
373
374 if ( verbose ) DenseMtx_writeForMatlab( mtxX, "mtxXX", matFile ) ;
375 if ( verbose ) fclose(matFile);
376
377 // Result->SolveTime().Time_First( ) ;
378 //
379 // Here we copy the values of B back from mtxX into xValues
380 //
381 if ( iam == 0 ) {
382 for ( int i = 0 ; i < numrows; i ++ )
383 {
384 DenseMtx_realEntry( mtxX, i, 0, &xValues[i] );
385 }
386 // DenseMtx_writeForMatlab( mtxX, "mtxX", matFile ) ;
387 }
388
389
390 if ( iam == 0 ) {
391 InpMtx_free(mtxA) ;
392 DenseMtx_free(mtxX) ;
393 DenseMtx_free(mtxY) ;
394 }
395 BridgeMPI_free(bridgeMPI) ;
396
397 Comm.Barrier();
398
399 return(0) ;
400}
401
static bool verbose
Definition Amesos.cpp:67
#define EPETRA_CHK_ERR(xxx)
Epetra_MultiVector * GetRHS() const
Definition SpoolesOO.h:70
int SetRHS(Epetra_MultiVector *B)
Definition SpoolesOO.cpp:99
Epetra_MultiVector * B_
Definition SpoolesOO.h:87
bool GetTrans() const
Definition SpoolesOO.h:72
Epetra_MultiVector * GetLHS() const
Definition SpoolesOO.h:68
virtual ~SpoolesOO(void)
Definition SpoolesOO.cpp:73
int Solve()
Epetra_Operator * UserOperator_
Definition SpoolesOO.h:82
double * x_
Definition SpoolesOO.h:93
int SetLHS(Epetra_MultiVector *X)
Definition SpoolesOO.cpp:90
int b_LDA_
Definition SpoolesOO.h:94
int SetSpoolesDefaults()
Epetra_RowMatrix * UserMatrix_
Definition SpoolesOO.h:83
bool inConstructor_
Definition SpoolesOO.h:96
Epetra_RowMatrix * GetUserMatrix() const
Definition SpoolesOO.h:66
int SetUserMatrix(Epetra_RowMatrix *UserMatrix)
Definition SpoolesOO.cpp:79
double * b_
Definition SpoolesOO.h:95
int x_LDA_
Definition SpoolesOO.h:92
Epetra_MultiVector * X_
Definition SpoolesOO.h:86