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 
00058 //
00059 //This is the Aztec implementation of LinSysCore.
00060 //
00061 namespace fei_trilinos {
00062 
00063 class Aztec_Map;
00064 class Aztec_BlockMap;
00065 class Aztec_LSVector;
00066 class AztecDMSR_Matrix;
00067 class AztecDVBR_Matrix;
00068 
00069 class Aztec_LinSysCore: public LinearSystemCore {
00070  public:
00071    Aztec_LinSysCore(MPI_Comm comm);
00072    virtual ~Aztec_LinSysCore();
00073 
00074    //for creating another instance of LinearSystemCore without knowing
00075    //the run-time type of 'this' one.
00076    LinearSystemCore* clone();
00077 
00078    //int parameters:
00079    //for setting generic argc/argv style parameters.
00080 
00081    int parameters(int numParams, const char*const * params);
00082 
00083    int setLookup(Lookup& lookup);
00084 
00085    int setGlobalOffsets(int len, int* nodeOffsets, int* eqnOffsets,
00086                                  int* blkEqnOffsets);
00087 
00088    int setConnectivities(GlobalID elemBlock,
00089                           int numElements,
00090                           int numNodesPerElem,
00091                           const GlobalID* elemIDs,
00092                           const int* const* connNodes) ;
00093 
00094    int setStiffnessMatrices(GlobalID,
00095           int,
00096           const GlobalID*,
00097           const double *const *const *,
00098           int,
00099           const int *const *)
00100      { return(0); }
00101 
00102    int setLoadVectors(GlobalID,
00103           int,
00104           const GlobalID*,
00105           const double *const *,
00106           int,
00107           const int *const *)
00108      { return(0); }
00109 
00110    int setMatrixStructure(int** ptColIndices,
00111                            int* ptRowLengths,
00112                            int** blkColIndices,
00113                            int* blkRowLengths,
00114                            int* ptRowsPerBlkRow) ;
00115 
00116    int setMultCREqns(int,
00117          int, int,
00118          int**, int**,
00119          int*,
00120          int*)
00121      { return(0); }
00122 
00123    int setPenCREqns(int,
00124          int, int,
00125          int**, int**,
00126          int*)
00127      { return(0); }
00128 
00129    //int resetMatrixAndVector:
00130    //don't destroy the structure of the matrix, but set the value 's'
00131    //throughout the matrix and vectors.
00132 
00133    int resetMatrixAndVector(double s);
00134    int resetMatrix(double s);
00135    int resetRHSVector(double s);
00136 
00137    //int sumIntoSystemMatrix:
00138    //this is the primary assembly function. The coefficients 'values'
00139    //are to be accumumlated into (added to any values already in place)
00140    //global (0-based) equations in 'ptRows' of the matrix.
00141 
00142    int sumIntoSystemMatrix(int numPtRows, const int* ptRows,
00143                             int numPtCols, const int* ptCols,
00144                             int numBlkRows, const int* blkRows,
00145                             int numBlkCols, const int* blkCols,
00146                             const double* const* values);
00147 
00148    int sumIntoSystemMatrix(int numPtRows, const int* ptRows,
00149                             int numPtCols, const int* ptCols,
00150                             const double* const* values);
00151    int putIntoSystemMatrix(int numPtRows, const int* ptRows,
00152                             int numPtCols, const int* ptCols,
00153                             const double* const* values);
00154 
00155    int getMatrixRowLength(int row, int& length);
00156 
00157    int getMatrixRow(int row, double* coefs, int* indices,
00158                             int len, int& rowLength);
00159 
00160    //int sumIntoRHSVector:
00161    //this is the rhs vector equivalent to sumIntoSystemMatrix above.
00162 
00163    int sumIntoRHSVector(int num,
00164                           const double* values,
00165                           const int* indices);
00166    int putIntoRHSVector(int num,
00167                           const double* values,
00168                           const int* indices);
00169    int getFromRHSVector(int num,
00170                           double* values,
00171                           const int* indices);
00172 
00173    //int matrixLoadComplete:
00174    //do any internal synchronization/communication.
00175 
00176    int matrixLoadComplete();
00177    
00178    //functions for enforcing boundary conditions.
00179    int enforceEssentialBC(int* globalEqn,
00180                            double* alpha,
00181                            double* gamma, int len);
00182 
00183    int enforceBlkEssentialBC(int* blkEqn, int* blkOffset,
00184                               double* alpha, double* gamma,
00185                               int len);
00186 
00187    int enforceRemoteEssBCs(int numEqns, int* globalEqns,
00188                             int** colIndices, int* colIndLen,
00189                             double** coefs);
00190 
00191    int enforceBlkRemoteEssBCs(int numEqns, int* blkEqns,
00192                                int** blkColInds, int** blkColOffsets,
00193                                int* blkColLens, double** remEssBCCoefs);
00194 
00195    //functions for getting/setting matrix or vector pointers.
00196 
00197    //getMatrixPtr:
00198    //obtain a pointer to the 'A' matrix. This should be considered a
00199    //constant pointer -- i.e., this class remains responsible for the
00200    //matrix (e.g., de-allocation upon destruction). 
00201    int getMatrixPtr(Data& data);
00202 
00203    //copyInMatrix:
00204    //replaces the internal matrix with a copy of the input argument, scaled
00205    //by the coefficient 'scalar'.
00206 
00207    int copyInMatrix(double scalar, const Data& data);
00208 
00209    //copyOutMatrix:
00210    //passes out a copy of the internal matrix, scaled by the coefficient
00211    //'scalar'.
00212 
00213    int copyOutMatrix(double scalar, Data& data);
00214 
00215    //sumInMatrix:
00216    //accumulate (sum) a copy of the input argument into the internal
00217    //matrix, scaling the input by the coefficient 'scalar'.
00218 
00219    int sumInMatrix(double scalar, const Data& data);
00220 
00221    //get/setRHSVectorPtr:
00222    //the same semantics apply here as for the matrixPtr functions above.
00223 
00224    int getRHSVectorPtr(Data& data);
00225 
00226    //copyInRHSVector/copyOutRHSVector/sumInRHSVector:
00227    //the same semantics apply here as for the matrix functions above.
00228 
00229    int copyInRHSVector(double scalar, const Data& data);
00230    int copyOutRHSVector(double scalar, Data& data);
00231    int sumInRHSVector(double scalar, const Data& data);
00232 
00233    //destroyMatrixData/destroyVectorData:
00234    //Utility function for destroying the matrix (or vector) in Data
00235 
00236    int destroyMatrixData(Data& data);
00237    int destroyVectorData(Data& data);
00238 
00239    //functions for managing multiple rhs vectors
00240    int setNumRHSVectors(int numRHSs, const int* rhsIDs);
00241 
00242    //int setRHSID:
00243    //set the 'current' rhs context, assuming multiple rhs vectors.
00244    int setRHSID(int rhsID);
00245 
00246    //int putInitialGuess:
00247    //function for setting (a subset of) the initial-guess
00248    //solution values (i.e., in the 'x' vector).
00249 
00250    int putInitialGuess(const int* eqnNumbers, const double* values,
00251                         int len);
00252 
00253    //function for getting all of the answers ('x' vector).
00254    int getSolution(double* answers, int len);
00255 
00256    //function for getting the (single) entry at equation
00257    //number 'eqnNumber'.
00258    int getSolnEntry(int eqnNumber, double& answer);
00259 
00260    //function for obtaining the residual vector (using current solution or
00261    //initial guess)
00262    int formResidual(double* values, int len);
00263 
00264    //function for launching the linear solver
00265    int launchSolver(int& solveStatus, int& iterations);
00266 
00267    int putNodalFieldData(int, int, int*, int, const double*)
00268      { return(0); }
00269 
00270    int writeSystem(const char* name);
00271 
00272  private:        //functions
00273 
00274    int createMiscStuff();
00275 
00276    int allocateMatrix(int** ptColIndices, int* ptRowLengths,
00277                        int** blkColIndices, int* blkRowLengths,
00278                        int* ptRowsPerBlkRow);
00279 
00280    int VBRmatPlusScaledMat(AztecDVBR_Matrix* A, double scalar,
00281                             AztecDVBR_Matrix* source);
00282 
00283    int MSRmatPlusScaledMat(AztecDMSR_Matrix* A, double scalar,
00284                             AztecDMSR_Matrix* source);
00285 
00286    int createBlockMatrix(int** blkColIndices,
00287                           int* blkRowLengths,
00288                           int* ptRowsPerBlkRow);
00289 
00290    int sumIntoBlockRow(int numBlkRows, const int* blkRows,
00291                         int numBlkCols, const int* blkCols,
00292                         const double* const* values,
00293                         int numPtCols,
00294            bool overwriteInsteadOfAccumulate);
00295 
00296    int copyBlockRow(int i, const int* blkRows,
00297                      int numBlkCols, const int* blkCols,
00298                      const double* const* values,
00299                      double* coefs);
00300 
00301    int modifyRHSforBCs();
00302 
00303    int explicitlySetDirichletBCs();
00304 
00305    int blockRowToPointRow(int blkRow);
00306 
00307    int getBlockRow(int blkRow, double*& val, int& valLen,
00308                    int*& blkColInds, int& blkColIndLen,
00309                    int& numNzBlks, int& numNNZ);
00310 
00311    int getBlkEqnsAndOffsets(int* ptEqns, int* blkEqns, int* blkOffsets,
00312                              int numEqns);
00313 
00314    int getBlockSize(int blkInd);
00315 
00316    int sumIntoPointRow(int numPtRows, const int* ptRows,
00317            int numPtCols, const int* ptColIndices,
00318            const double* const* values,
00319            bool overwriteInsteadOfAccumulate);
00320 
00321    int sumPointIntoBlockRow(int blkRow, int rowOffset,
00322           int blkCol, int colOffset, double value);
00323 
00324    int setMatrixType(const char* name);
00325    int selectSolver(const char* name);
00326    int selectPreconditioner(const char* name);
00327    void setSubdomainSolve(const char* name);
00328    void setScalingOption(const char* param);
00329    void setConvTest(const char* param);
00330    void setPreCalc(const char* param);
00331    void setTypeOverlap(const char* param);
00332    void setOverlap(const char* param);
00333    void setOrthog(const char* param);
00334    void setAuxVec(const char* param);
00335    void setAZ_output(const char* param);
00336 
00337    void recordUserParams();
00338 
00339    void checkForParam(const char* paramName, int numParams_,
00340                       char** paramStrings,
00341                       double& param);
00342 
00343    void recordUserOptions();
00344 
00345    void checkForOption(const char* paramName, int numParams_,
00346                        char** paramStrings,
00347                        int& param);
00348 
00349    int blkRowEssBCMod(int blkEqn, int blkOffset, double* val, int* blkCols,
00350                         int numCols, int numPtNNZ, double alpha, double gamma);
00351 
00352    int blkColEssBCMod(int blkRow, int blkEqn, int blkOffset, double* val,
00353                       int* blkCols,
00354                       int numCols, int numPtNNZ, double alpha, double gamma);
00355 
00356    void setDebugOutput(const char* path, const char* name);
00357 
00358    void debugOutput(const char* msg) const;
00359 
00360    int writeA(const char* name);
00361    int writeVec(Aztec_LSVector* v, const char* name);
00362 
00363    int messageAbort(const char* msg) const;
00364 
00365  private:            //variables
00366 
00367    MPI_Comm comm_;
00368 
00369    Lookup* lookup_;
00370    bool haveLookup_;
00371 
00372    int numProcs_;
00373    int thisProc_;
00374    int masterProc_;
00375 
00376    int* update_;
00377    fei::SharedPtr<Aztec_Map> map_;
00378    AztecDMSR_Matrix *A_;
00379    AztecDMSR_Matrix *A_ptr_;
00380    Aztec_LSVector *x_, **b_, *bc_;
00381    int* essBCindices_;
00382    int numEssBCs_;
00383    bool bcsLoaded_;
00384    bool explicitDirichletBCs_;
00385    bool BCenforcement_no_column_mod_;
00386    Aztec_LSVector *b_ptr_;
00387    bool matrixAllocated_;
00388    bool vectorsAllocated_;
00389    bool blkMatrixAllocated_;
00390    bool matrixLoaded_;
00391    bool rhsLoaded_;
00392    bool needNewPreconditioner_;
00393 
00394    bool tooLateToChooseBlock_;
00395    bool blockMatrix_;
00396    fei::SharedPtr<Aztec_BlockMap> blkMap_;
00397    AztecDVBR_Matrix *blkA_;
00398    AztecDVBR_Matrix *blkA_ptr_;
00399    int* blkUpdate_;
00400 
00401    AZ_MATRIX *azA_;
00402    AZ_PRECOND *azP_;
00403    bool precondCreated_;
00404    AZ_SCALING *azS_;
00405    bool scalingCreated_;
00406 
00407    int *aztec_options_;
00408    double *aztec_params_;
00409    double *aztec_status_;
00410 
00411    double* tmp_x_;
00412    bool tmp_x_touched_;
00413    double** tmp_b_;
00414    double* tmp_bc_;
00415    bool tmp_b_allocated_;
00416 
00417    bool ML_Vanek_;
00418    int numLevels_;
00419 
00420    int* rhsIDs_;
00421    int numRHSs_;
00422 
00423    int currentRHS_;
00424 
00425    int numGlobalEqns_;
00426    int localOffset_;
00427    int numLocalEqns_;
00428 
00429    int numGlobalEqnBlks_;
00430    int localBlkOffset_;
00431    int numLocalEqnBlks_;
00432    int* localBlockSizes_;
00433 
00434    int numNonzeroBlocks_;
00435 
00436    int outputLevel_;
00437    int numParams_;
00438    char** paramStrings_;
00439 
00440    std::string name_;
00441    int debugOutput_;
00442    int debugFileCounter_;
00443    char* debugPath_;
00444    char* debugFileName_;
00445    FILE* debugFile_;
00446 
00447    std::map<std::string,unsigned>& named_solve_counter_;
00448 };
00449 
00450 }//namespace fei_trilinos
00451 
00452 #endif
00453 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends