FEI Version of the Day
Loading...
Searching...
No Matches
fei_Aztec_LinSysCore.hpp
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43#ifndef _fei_Aztec_LinSysCore_hpp_
44#define _fei_Aztec_LinSysCore_hpp_
45
46
47#include <fei_macros.hpp>
48#include <fei_defs.h>
49#include <fei_Data.hpp>
50#include <fei_LinearSystemCore.hpp>
51#include <fei_SharedPtr.hpp>
52
53#include <string>
54#include <map>
55
56#include <az_aztec.h>
57#include <fei_AztecDMSR_Matrix.hpp>
58
59//
60//This is the Aztec implementation of LinSysCore.
61//
62namespace fei_trilinos {
63
64class Aztec_Map;
65class Aztec_BlockMap;
66class Aztec_LSVector;
67class AztecDMSR_Matrix;
68class AztecDVBR_Matrix;
69
70class Aztec_LinSysCore: public LinearSystemCore {
71 public:
72 Aztec_LinSysCore(MPI_Comm comm);
73 virtual ~Aztec_LinSysCore();
74
75 //for creating another instance of LinearSystemCore without knowing
76 //the run-time type of 'this' one.
77 LinearSystemCore* clone();
78
79 //int parameters:
80 //for setting generic argc/argv style parameters.
81
82 int parameters(int numParams, const char*const * params);
83
84 int setLookup(Lookup& lookup);
85
86 int setGlobalOffsets(int len, int* nodeOffsets, int* eqnOffsets,
87 int* blkEqnOffsets);
88
89 int setConnectivities(GlobalID elemBlock,
90 int numElements,
91 int numNodesPerElem,
92 const GlobalID* elemIDs,
93 const int* const* connNodes) ;
94
95 int setStiffnessMatrices(GlobalID,
96 int,
97 const GlobalID*,
98 const double *const *const *,
99 int,
100 const int *const *)
101 { return(0); }
102
103 int setLoadVectors(GlobalID,
104 int,
105 const GlobalID*,
106 const double *const *,
107 int,
108 const int *const *)
109 { return(0); }
110
111 int setMatrixStructure(int** ptColIndices,
112 int* ptRowLengths,
113 int** blkColIndices,
114 int* blkRowLengths,
115 int* ptRowsPerBlkRow) ;
116
117 int setMultCREqns(int,
118 int, int,
119 int**, int**,
120 int*,
121 int*)
122 { return(0); }
123
124 int setPenCREqns(int,
125 int, int,
126 int**, int**,
127 int*)
128 { return(0); }
129
130 //int resetMatrixAndVector:
131 //don't destroy the structure of the matrix, but set the value 's'
132 //throughout the matrix and vectors.
133
134 int resetMatrixAndVector(double s);
135 int resetMatrix(double s);
136 int resetRHSVector(double s);
137
138 //int sumIntoSystemMatrix:
139 //this is the primary assembly function. The coefficients 'values'
140 //are to be accumumlated into (added to any values already in place)
141 //global (0-based) equations in 'ptRows' of the matrix.
142
143 int sumIntoSystemMatrix(int numPtRows, const int* ptRows,
144 int numPtCols, const int* ptCols,
145 int numBlkRows, const int* blkRows,
146 int numBlkCols, const int* blkCols,
147 const double* const* values);
148
149 int sumIntoSystemMatrix(int numPtRows, const int* ptRows,
150 int numPtCols, const int* ptCols,
151 const double* const* values);
152 int putIntoSystemMatrix(int numPtRows, const int* ptRows,
153 int numPtCols, const int* ptCols,
154 const double* const* values);
155
156 int getMatrixRowLength(int row, int& length);
157
158 int getMatrixRow(int row, double* coefs, int* indices,
159 int len, int& rowLength);
160
161 //int sumIntoRHSVector:
162 //this is the rhs vector equivalent to sumIntoSystemMatrix above.
163
164 int sumIntoRHSVector(int num,
165 const double* values,
166 const int* indices);
167 int putIntoRHSVector(int num,
168 const double* values,
169 const int* indices);
170 int getFromRHSVector(int num,
171 double* values,
172 const int* indices);
173
174 //int matrixLoadComplete:
175 //do any internal synchronization/communication.
176
177 int matrixLoadComplete();
178
179 //functions for enforcing boundary conditions.
180 int enforceEssentialBC(int* globalEqn,
181 double* alpha,
182 double* gamma, int len);
183
184 int enforceBlkEssentialBC(int* blkEqn, int* blkOffset,
185 double* alpha, double* gamma,
186 int len);
187
188 int enforceRemoteEssBCs(int numEqns, int* globalEqns,
189 int** colIndices, int* colIndLen,
190 double** coefs);
191
192 int enforceBlkRemoteEssBCs(int numEqns, int* blkEqns,
193 int** blkColInds, int** blkColOffsets,
194 int* blkColLens, double** remEssBCCoefs);
195
196 //functions for getting/setting matrix or vector pointers.
197
198 //getMatrixPtr:
199 //obtain a pointer to the 'A' matrix. This should be considered a
200 //constant pointer -- i.e., this class remains responsible for the
201 //matrix (e.g., de-allocation upon destruction).
202 int getMatrixPtr(Data& data);
203
204 //copyInMatrix:
205 //replaces the internal matrix with a copy of the input argument, scaled
206 //by the coefficient 'scalar'.
207
208 int copyInMatrix(double scalar, const Data& data);
209
210 //copyOutMatrix:
211 //passes out a copy of the internal matrix, scaled by the coefficient
212 //'scalar'.
213
214 int copyOutMatrix(double scalar, Data& data);
215
216 //sumInMatrix:
217 //accumulate (sum) a copy of the input argument into the internal
218 //matrix, scaling the input by the coefficient 'scalar'.
219
220 int sumInMatrix(double scalar, const Data& data);
221
222 //get/setRHSVectorPtr:
223 //the same semantics apply here as for the matrixPtr functions above.
224
225 int getRHSVectorPtr(Data& data);
226
227 //copyInRHSVector/copyOutRHSVector/sumInRHSVector:
228 //the same semantics apply here as for the matrix functions above.
229
230 int copyInRHSVector(double scalar, const Data& data);
231 int copyOutRHSVector(double scalar, Data& data);
232 int sumInRHSVector(double scalar, const Data& data);
233
234 //destroyMatrixData/destroyVectorData:
235 //Utility function for destroying the matrix (or vector) in Data
236
237 int destroyMatrixData(Data& data);
238 int destroyVectorData(Data& data);
239
240 //functions for managing multiple rhs vectors
241 int setNumRHSVectors(int numRHSs, const int* rhsIDs);
242
243 //int setRHSID:
244 //set the 'current' rhs context, assuming multiple rhs vectors.
245 int setRHSID(int rhsID);
246
247 //int putInitialGuess:
248 //function for setting (a subset of) the initial-guess
249 //solution values (i.e., in the 'x' vector).
250
251 int putInitialGuess(const int* eqnNumbers, const double* values,
252 int len);
253
254 //function for getting all of the answers ('x' vector).
255 int getSolution(double* answers, int len);
256
257 //function for getting the (single) entry at equation
258 //number 'eqnNumber'.
259 int getSolnEntry(int eqnNumber, double& answer);
260
261 //function for obtaining the residual vector (using current solution or
262 //initial guess)
263 int formResidual(double* values, int len);
264
265 //function for launching the linear solver
266 int launchSolver(int& solveStatus, int& iterations);
267
268 int putNodalFieldData(int, int, int*, int, const double*)
269 { return(0); }
270
271 int writeSystem(const char* name);
272
273 double* getMatrixBeginPointer() { return A_ptr_->getBeginPointer(); }
274
275 int getMatrixOffset(int row, int col) { return A_ptr_->getOffset(row,col); }
276
277 private: //functions
278
279 int createMiscStuff();
280
281 int allocateMatrix(int** ptColIndices, int* ptRowLengths,
282 int** blkColIndices, int* blkRowLengths,
283 int* ptRowsPerBlkRow);
284
285 int VBRmatPlusScaledMat(AztecDVBR_Matrix* A, double scalar,
286 AztecDVBR_Matrix* source);
287
288 int MSRmatPlusScaledMat(AztecDMSR_Matrix* A, double scalar,
289 AztecDMSR_Matrix* source);
290
291 int createBlockMatrix(int** blkColIndices,
292 int* blkRowLengths,
293 int* ptRowsPerBlkRow);
294
295 int sumIntoBlockRow(int numBlkRows, const int* blkRows,
296 int numBlkCols, const int* blkCols,
297 const double* const* values,
298 int numPtCols,
299 bool overwriteInsteadOfAccumulate);
300
301 int copyBlockRow(int i, const int* blkRows,
302 int numBlkCols, const int* blkCols,
303 const double* const* values,
304 double* coefs);
305
306 int modifyRHSforBCs();
307
308 int explicitlySetDirichletBCs();
309
310 int blockRowToPointRow(int blkRow);
311
312 int getBlockRow(int blkRow, double*& val, int& valLen,
313 int*& blkColInds, int& blkColIndLen,
314 int& numNzBlks, int& numNNZ);
315
316 int getBlkEqnsAndOffsets(int* ptEqns, int* blkEqns, int* blkOffsets,
317 int numEqns);
318
319 int getBlockSize(int blkInd);
320
321 int sumIntoPointRow(int numPtRows, const int* ptRows,
322 int numPtCols, const int* ptColIndices,
323 const double* const* values,
324 bool overwriteInsteadOfAccumulate);
325
326 int sumPointIntoBlockRow(int blkRow, int rowOffset,
327 int blkCol, int colOffset, double value);
328
329 int setMatrixType(const char* name);
330 int selectSolver(const char* name);
331 int selectPreconditioner(const char* name);
332 void setSubdomainSolve(const char* name);
333 void setScalingOption(const char* param);
334 void setConvTest(const char* param);
335 void setPreCalc(const char* param);
336 void setTypeOverlap(const char* param);
337 void setOverlap(const char* param);
338 void setOrthog(const char* param);
339 void setAuxVec(const char* param);
340 void setAZ_output(const char* param);
341
342 void recordUserParams();
343
344 void checkForParam(const char* paramName, int numParams_,
345 char** paramStrings,
346 double& param);
347
348 void recordUserOptions();
349
350 void checkForOption(const char* paramName, int numParams_,
351 char** paramStrings,
352 int& param);
353
354 int blkRowEssBCMod(int blkEqn, int blkOffset, double* val, int* blkCols,
355 int numCols, int numPtNNZ, double alpha, double gamma);
356
357 int blkColEssBCMod(int blkRow, int blkEqn, int blkOffset, double* val,
358 int* blkCols,
359 int numCols, int numPtNNZ, double alpha, double gamma);
360
361 void setDebugOutput(const char* path, const char* name);
362
363 void debugOutput(const char* msg) const;
364
365 int writeA(const char* name);
366 int writeVec(Aztec_LSVector* v, const char* name);
367
368 int messageAbort(const char* msg) const;
369
370 private: //variables
371
372 MPI_Comm comm_;
373
374 Lookup* lookup_;
375 bool haveLookup_;
376
377 int numProcs_;
378 int thisProc_;
379 int masterProc_;
380
381 int* update_;
383 AztecDMSR_Matrix *A_;
384 AztecDMSR_Matrix *A_ptr_;
385 Aztec_LSVector *x_, **b_, *bc_;
386 int* essBCindices_;
387 int numEssBCs_;
388 bool bcsLoaded_;
389 bool explicitDirichletBCs_;
390 bool BCenforcement_no_column_mod_;
391 Aztec_LSVector *b_ptr_;
392 bool matrixAllocated_;
393 bool vectorsAllocated_;
394 bool blkMatrixAllocated_;
395 bool matrixLoaded_;
396 bool rhsLoaded_;
397 bool needNewPreconditioner_;
398
399 bool tooLateToChooseBlock_;
400 bool blockMatrix_;
402 AztecDVBR_Matrix *blkA_;
403 AztecDVBR_Matrix *blkA_ptr_;
404 int* blkUpdate_;
405
406 AZ_MATRIX *azA_;
407 AZ_PRECOND *azP_;
408 bool precondCreated_;
409 AZ_SCALING *azS_;
410 bool scalingCreated_;
411
412 int *aztec_options_;
413 double *aztec_params_;
414 double *aztec_status_;
415
416 double* tmp_x_;
417 bool tmp_x_touched_;
418 double** tmp_b_;
419 double* tmp_bc_;
420 bool tmp_b_allocated_;
421
422 bool ML_Vanek_;
423 int numLevels_;
424
425 int* rhsIDs_;
426 int numRHSs_;
427
428 int currentRHS_;
429
430 int numGlobalEqns_;
431 int localOffset_;
432 int numLocalEqns_;
433
434 int numGlobalEqnBlks_;
435 int localBlkOffset_;
436 int numLocalEqnBlks_;
437 int* localBlockSizes_;
438
439 int numNonzeroBlocks_;
440
441 int outputLevel_;
442 int numParams_;
443 char** paramStrings_;
444
445 std::string name_;
446 int debugOutput_;
447 int debugFileCounter_;
448 char* debugPath_;
449 char* debugFileName_;
450 FILE* debugFile_;
451
452 std::map<std::string,unsigned>& named_solve_counter_;
453};
454
455}//namespace fei_trilinos
456
457#endif
458