FEI Version of the Day
fei_Aztec_LinSysCore.hpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 #ifndef _fei_Aztec_LinSysCore_hpp_
00044 #define _fei_Aztec_LinSysCore_hpp_
00045 
00046 
00047 #include <fei_macros.hpp>
00048 #include <fei_defs.h>
00049 #include <fei_Data.hpp>
00050 #include <fei_LinearSystemCore.hpp>
00051 #include <fei_SharedPtr.hpp>
00052 
00053 #include <string>
00054 #include <map>
00055 
00056 #include <az_aztec.h>
00057 #include <fei_AztecDMSR_Matrix.hpp>
00058 
00059 //
00060 //This is the Aztec implementation of LinSysCore.
00061 //
00062 namespace fei_trilinos {
00063 
00064 class Aztec_Map;
00065 class Aztec_BlockMap;
00066 class Aztec_LSVector;
00067 class AztecDMSR_Matrix;
00068 class AztecDVBR_Matrix;
00069 
00070 class Aztec_LinSysCore: public LinearSystemCore {
00071  public:
00072    Aztec_LinSysCore(MPI_Comm comm);
00073    virtual ~Aztec_LinSysCore();
00074 
00075    //for creating another instance of LinearSystemCore without knowing
00076    //the run-time type of 'this' one.
00077    LinearSystemCore* clone();
00078 
00079    //int parameters:
00080    //for setting generic argc/argv style parameters.
00081 
00082    int parameters(int numParams, const char*const * params);
00083 
00084    int setLookup(Lookup& lookup);
00085 
00086    int setGlobalOffsets(int len, int* nodeOffsets, int* eqnOffsets,
00087                                  int* blkEqnOffsets);
00088 
00089    int setConnectivities(GlobalID elemBlock,
00090                           int numElements,
00091                           int numNodesPerElem,
00092                           const GlobalID* elemIDs,
00093                           const int* const* connNodes) ;
00094 
00095    int setStiffnessMatrices(GlobalID,
00096           int,
00097           const GlobalID*,
00098           const double *const *const *,
00099           int,
00100           const int *const *)
00101      { return(0); }
00102 
00103    int setLoadVectors(GlobalID,
00104           int,
00105           const GlobalID*,
00106           const double *const *,
00107           int,
00108           const int *const *)
00109      { return(0); }
00110 
00111    int setMatrixStructure(int** ptColIndices,
00112                            int* ptRowLengths,
00113                            int** blkColIndices,
00114                            int* blkRowLengths,
00115                            int* ptRowsPerBlkRow) ;
00116 
00117    int setMultCREqns(int,
00118          int, int,
00119          int**, int**,
00120          int*,
00121          int*)
00122      { return(0); }
00123 
00124    int setPenCREqns(int,
00125          int, int,
00126          int**, int**,
00127          int*)
00128      { return(0); }
00129 
00130    //int resetMatrixAndVector:
00131    //don't destroy the structure of the matrix, but set the value 's'
00132    //throughout the matrix and vectors.
00133 
00134    int resetMatrixAndVector(double s);
00135    int resetMatrix(double s);
00136    int resetRHSVector(double s);
00137 
00138    //int sumIntoSystemMatrix:
00139    //this is the primary assembly function. The coefficients 'values'
00140    //are to be accumumlated into (added to any values already in place)
00141    //global (0-based) equations in 'ptRows' of the matrix.
00142 
00143    int sumIntoSystemMatrix(int numPtRows, const int* ptRows,
00144                             int numPtCols, const int* ptCols,
00145                             int numBlkRows, const int* blkRows,
00146                             int numBlkCols, const int* blkCols,
00147                             const double* const* values);
00148 
00149    int sumIntoSystemMatrix(int numPtRows, const int* ptRows,
00150                             int numPtCols, const int* ptCols,
00151                             const double* const* values);
00152    int putIntoSystemMatrix(int numPtRows, const int* ptRows,
00153                             int numPtCols, const int* ptCols,
00154                             const double* const* values);
00155 
00156    int getMatrixRowLength(int row, int& length);
00157 
00158    int getMatrixRow(int row, double* coefs, int* indices,
00159                             int len, int& rowLength);
00160 
00161    //int sumIntoRHSVector:
00162    //this is the rhs vector equivalent to sumIntoSystemMatrix above.
00163 
00164    int sumIntoRHSVector(int num,
00165                           const double* values,
00166                           const int* indices);
00167    int putIntoRHSVector(int num,
00168                           const double* values,
00169                           const int* indices);
00170    int getFromRHSVector(int num,
00171                           double* values,
00172                           const int* indices);
00173 
00174    //int matrixLoadComplete:
00175    //do any internal synchronization/communication.
00176 
00177    int matrixLoadComplete();
00178    
00179    //functions for enforcing boundary conditions.
00180    int enforceEssentialBC(int* globalEqn,
00181                            double* alpha,
00182                            double* gamma, int len);
00183 
00184    int enforceBlkEssentialBC(int* blkEqn, int* blkOffset,
00185                               double* alpha, double* gamma,
00186                               int len);
00187 
00188    int enforceRemoteEssBCs(int numEqns, int* globalEqns,
00189                             int** colIndices, int* colIndLen,
00190                             double** coefs);
00191 
00192    int enforceBlkRemoteEssBCs(int numEqns, int* blkEqns,
00193                                int** blkColInds, int** blkColOffsets,
00194                                int* blkColLens, double** remEssBCCoefs);
00195 
00196    //functions for getting/setting matrix or vector pointers.
00197 
00198    //getMatrixPtr:
00199    //obtain a pointer to the 'A' matrix. This should be considered a
00200    //constant pointer -- i.e., this class remains responsible for the
00201    //matrix (e.g., de-allocation upon destruction). 
00202    int getMatrixPtr(Data& data);
00203 
00204    //copyInMatrix:
00205    //replaces the internal matrix with a copy of the input argument, scaled
00206    //by the coefficient 'scalar'.
00207 
00208    int copyInMatrix(double scalar, const Data& data);
00209 
00210    //copyOutMatrix:
00211    //passes out a copy of the internal matrix, scaled by the coefficient
00212    //'scalar'.
00213 
00214    int copyOutMatrix(double scalar, Data& data);
00215 
00216    //sumInMatrix:
00217    //accumulate (sum) a copy of the input argument into the internal
00218    //matrix, scaling the input by the coefficient 'scalar'.
00219 
00220    int sumInMatrix(double scalar, const Data& data);
00221 
00222    //get/setRHSVectorPtr:
00223    //the same semantics apply here as for the matrixPtr functions above.
00224 
00225    int getRHSVectorPtr(Data& data);
00226 
00227    //copyInRHSVector/copyOutRHSVector/sumInRHSVector:
00228    //the same semantics apply here as for the matrix functions above.
00229 
00230    int copyInRHSVector(double scalar, const Data& data);
00231    int copyOutRHSVector(double scalar, Data& data);
00232    int sumInRHSVector(double scalar, const Data& data);
00233 
00234    //destroyMatrixData/destroyVectorData:
00235    //Utility function for destroying the matrix (or vector) in Data
00236 
00237    int destroyMatrixData(Data& data);
00238    int destroyVectorData(Data& data);
00239 
00240    //functions for managing multiple rhs vectors
00241    int setNumRHSVectors(int numRHSs, const int* rhsIDs);
00242 
00243    //int setRHSID:
00244    //set the 'current' rhs context, assuming multiple rhs vectors.
00245    int setRHSID(int rhsID);
00246 
00247    //int putInitialGuess:
00248    //function for setting (a subset of) the initial-guess
00249    //solution values (i.e., in the 'x' vector).
00250 
00251    int putInitialGuess(const int* eqnNumbers, const double* values,
00252                         int len);
00253 
00254    //function for getting all of the answers ('x' vector).
00255    int getSolution(double* answers, int len);
00256 
00257    //function for getting the (single) entry at equation
00258    //number 'eqnNumber'.
00259    int getSolnEntry(int eqnNumber, double& answer);
00260 
00261    //function for obtaining the residual vector (using current solution or
00262    //initial guess)
00263    int formResidual(double* values, int len);
00264 
00265    //function for launching the linear solver
00266    int launchSolver(int& solveStatus, int& iterations);
00267 
00268    int putNodalFieldData(int, int, int*, int, const double*)
00269      { return(0); }
00270 
00271    int writeSystem(const char* name);
00272 
00273    double* getMatrixBeginPointer() { return A_ptr_->getBeginPointer(); }
00274 
00275    int getMatrixOffset(int row, int col) { return A_ptr_->getOffset(row,col); }
00276 
00277  private:        //functions
00278 
00279    int createMiscStuff();
00280 
00281    int allocateMatrix(int** ptColIndices, int* ptRowLengths,
00282                        int** blkColIndices, int* blkRowLengths,
00283                        int* ptRowsPerBlkRow);
00284 
00285    int VBRmatPlusScaledMat(AztecDVBR_Matrix* A, double scalar,
00286                             AztecDVBR_Matrix* source);
00287 
00288    int MSRmatPlusScaledMat(AztecDMSR_Matrix* A, double scalar,
00289                             AztecDMSR_Matrix* source);
00290 
00291    int createBlockMatrix(int** blkColIndices,
00292                           int* blkRowLengths,
00293                           int* ptRowsPerBlkRow);
00294 
00295    int sumIntoBlockRow(int numBlkRows, const int* blkRows,
00296                         int numBlkCols, const int* blkCols,
00297                         const double* const* values,
00298                         int numPtCols,
00299            bool overwriteInsteadOfAccumulate);
00300 
00301    int copyBlockRow(int i, const int* blkRows,
00302                      int numBlkCols, const int* blkCols,
00303                      const double* const* values,
00304                      double* coefs);
00305 
00306    int modifyRHSforBCs();
00307 
00308    int explicitlySetDirichletBCs();
00309 
00310    int blockRowToPointRow(int blkRow);
00311 
00312    int getBlockRow(int blkRow, double*& val, int& valLen,
00313                    int*& blkColInds, int& blkColIndLen,
00314                    int& numNzBlks, int& numNNZ);
00315 
00316    int getBlkEqnsAndOffsets(int* ptEqns, int* blkEqns, int* blkOffsets,
00317                              int numEqns);
00318 
00319    int getBlockSize(int blkInd);
00320 
00321    int sumIntoPointRow(int numPtRows, const int* ptRows,
00322            int numPtCols, const int* ptColIndices,
00323            const double* const* values,
00324            bool overwriteInsteadOfAccumulate);
00325 
00326    int sumPointIntoBlockRow(int blkRow, int rowOffset,
00327           int blkCol, int colOffset, double value);
00328 
00329    int setMatrixType(const char* name);
00330    int selectSolver(const char* name);
00331    int selectPreconditioner(const char* name);
00332    void setSubdomainSolve(const char* name);
00333    void setScalingOption(const char* param);
00334    void setConvTest(const char* param);
00335    void setPreCalc(const char* param);
00336    void setTypeOverlap(const char* param);
00337    void setOverlap(const char* param);
00338    void setOrthog(const char* param);
00339    void setAuxVec(const char* param);
00340    void setAZ_output(const char* param);
00341 
00342    void recordUserParams();
00343 
00344    void checkForParam(const char* paramName, int numParams_,
00345                       char** paramStrings,
00346                       double& param);
00347 
00348    void recordUserOptions();
00349 
00350    void checkForOption(const char* paramName, int numParams_,
00351                        char** paramStrings,
00352                        int& param);
00353 
00354    int blkRowEssBCMod(int blkEqn, int blkOffset, double* val, int* blkCols,
00355                         int numCols, int numPtNNZ, double alpha, double gamma);
00356 
00357    int blkColEssBCMod(int blkRow, int blkEqn, int blkOffset, double* val,
00358                       int* blkCols,
00359                       int numCols, int numPtNNZ, double alpha, double gamma);
00360 
00361    void setDebugOutput(const char* path, const char* name);
00362 
00363    void debugOutput(const char* msg) const;
00364 
00365    int writeA(const char* name);
00366    int writeVec(Aztec_LSVector* v, const char* name);
00367 
00368    int messageAbort(const char* msg) const;
00369 
00370  private:            //variables
00371 
00372    MPI_Comm comm_;
00373 
00374    Lookup* lookup_;
00375    bool haveLookup_;
00376 
00377    int numProcs_;
00378    int thisProc_;
00379    int masterProc_;
00380 
00381    int* update_;
00382    fei::SharedPtr<Aztec_Map> map_;
00383    AztecDMSR_Matrix *A_;
00384    AztecDMSR_Matrix *A_ptr_;
00385    Aztec_LSVector *x_, **b_, *bc_;
00386    int* essBCindices_;
00387    int numEssBCs_;
00388    bool bcsLoaded_;
00389    bool explicitDirichletBCs_;
00390    bool BCenforcement_no_column_mod_;
00391    Aztec_LSVector *b_ptr_;
00392    bool matrixAllocated_;
00393    bool vectorsAllocated_;
00394    bool blkMatrixAllocated_;
00395    bool matrixLoaded_;
00396    bool rhsLoaded_;
00397    bool needNewPreconditioner_;
00398 
00399    bool tooLateToChooseBlock_;
00400    bool blockMatrix_;
00401    fei::SharedPtr<Aztec_BlockMap> blkMap_;
00402    AztecDVBR_Matrix *blkA_;
00403    AztecDVBR_Matrix *blkA_ptr_;
00404    int* blkUpdate_;
00405 
00406    AZ_MATRIX *azA_;
00407    AZ_PRECOND *azP_;
00408    bool precondCreated_;
00409    AZ_SCALING *azS_;
00410    bool scalingCreated_;
00411 
00412    int *aztec_options_;
00413    double *aztec_params_;
00414    double *aztec_status_;
00415 
00416    double* tmp_x_;
00417    bool tmp_x_touched_;
00418    double** tmp_b_;
00419    double* tmp_bc_;
00420    bool tmp_b_allocated_;
00421 
00422    bool ML_Vanek_;
00423    int numLevels_;
00424 
00425    int* rhsIDs_;
00426    int numRHSs_;
00427 
00428    int currentRHS_;
00429 
00430    int numGlobalEqns_;
00431    int localOffset_;
00432    int numLocalEqns_;
00433 
00434    int numGlobalEqnBlks_;
00435    int localBlkOffset_;
00436    int numLocalEqnBlks_;
00437    int* localBlockSizes_;
00438 
00439    int numNonzeroBlocks_;
00440 
00441    int outputLevel_;
00442    int numParams_;
00443    char** paramStrings_;
00444 
00445    std::string name_;
00446    int debugOutput_;
00447    int debugFileCounter_;
00448    char* debugPath_;
00449    char* debugFileName_;
00450    FILE* debugFile_;
00451 
00452    std::map<std::string,unsigned>& named_solve_counter_;
00453 };
00454 
00455 }//namespace fei_trilinos
00456 
00457 #endif
00458 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends