fei_Aztec_LinSysCore.hpp

/*
// @HEADER
// ************************************************************************
//             FEI: Finite Element Interface to Linear Solvers
//                  Copyright (2005) Sandia Corporation.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
// U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Alan Williams (william@sandia.gov) 
//
// ************************************************************************
// @HEADER
*/

#ifndef _fei_Aztec_LinSysCore_hpp_
#define _fei_Aztec_LinSysCore_hpp_


#include <fei_macros.hpp>
#include <fei_defs.h>
#include <fei_Data.hpp>
#include <fei_LinearSystemCore.hpp>
#include <fei_SharedPtr.hpp>

#include <string>
#include <map>

#include <az_aztec.h>
#include <fei_AztecDMSR_Matrix.hpp>

//
//This is the Aztec implementation of LinSysCore.
//
namespace fei_trilinos {

class Aztec_Map;
class Aztec_BlockMap;
class Aztec_LSVector;
class AztecDMSR_Matrix;
class AztecDVBR_Matrix;

class Aztec_LinSysCore: public LinearSystemCore {
 public:
   Aztec_LinSysCore(MPI_Comm comm);
   virtual ~Aztec_LinSysCore();

   //for creating another instance of LinearSystemCore without knowing
   //the run-time type of 'this' one.
   LinearSystemCore* clone();

   //int parameters:
   //for setting generic argc/argv style parameters.

   int parameters(int numParams, const char*const * params);

   int setLookup(Lookup& lookup);

   int setGlobalOffsets(int len, int* nodeOffsets, int* eqnOffsets,
                                 int* blkEqnOffsets);

   int setConnectivities(GlobalID elemBlock,
                          int numElements,
                          int numNodesPerElem,
                          const GlobalID* elemIDs,
                          const int* const* connNodes) ;

   int setStiffnessMatrices(GlobalID,
          int,
          const GlobalID*,
          const double *const *const *,
          int,
          const int *const *)
     { return(0); }

   int setLoadVectors(GlobalID,
          int,
          const GlobalID*,
          const double *const *,
          int,
          const int *const *)
     { return(0); }

   int setMatrixStructure(int** ptColIndices,
                           int* ptRowLengths,
                           int** blkColIndices,
                           int* blkRowLengths,
                           int* ptRowsPerBlkRow) ;

   int setMultCREqns(int,
         int, int,
         int**, int**,
         int*,
         int*)
     { return(0); }

   int setPenCREqns(int,
         int, int,
         int**, int**,
         int*)
     { return(0); }

   //int resetMatrixAndVector:
   //don't destroy the structure of the matrix, but set the value 's'
   //throughout the matrix and vectors.

   int resetMatrixAndVector(double s);
   int resetMatrix(double s);
   int resetRHSVector(double s);

   //int sumIntoSystemMatrix:
   //this is the primary assembly function. The coefficients 'values'
   //are to be accumumlated into (added to any values already in place)
   //global (0-based) equations in 'ptRows' of the matrix.

   int sumIntoSystemMatrix(int numPtRows, const int* ptRows,
                            int numPtCols, const int* ptCols,
                            int numBlkRows, const int* blkRows,
                            int numBlkCols, const int* blkCols,
                            const double* const* values);

   int sumIntoSystemMatrix(int numPtRows, const int* ptRows,
                            int numPtCols, const int* ptCols,
                            const double* const* values);
   int putIntoSystemMatrix(int numPtRows, const int* ptRows,
                            int numPtCols, const int* ptCols,
                            const double* const* values);

   int getMatrixRowLength(int row, int& length);

   int getMatrixRow(int row, double* coefs, int* indices,
                            int len, int& rowLength);

   //int sumIntoRHSVector:
   //this is the rhs vector equivalent to sumIntoSystemMatrix above.

   int sumIntoRHSVector(int num,
                          const double* values,
                          const int* indices);
   int putIntoRHSVector(int num,
                          const double* values,
                          const int* indices);
   int getFromRHSVector(int num,
                          double* values,
                          const int* indices);

   //int matrixLoadComplete:
   //do any internal synchronization/communication.

   int matrixLoadComplete();
   
   //functions for enforcing boundary conditions.
   int enforceEssentialBC(int* globalEqn,
                           double* alpha,
                           double* gamma, int len);

   int enforceBlkEssentialBC(int* blkEqn, int* blkOffset,
                              double* alpha, double* gamma,
                              int len);

   int enforceRemoteEssBCs(int numEqns, int* globalEqns,
                            int** colIndices, int* colIndLen,
                            double** coefs);

   int enforceBlkRemoteEssBCs(int numEqns, int* blkEqns,
                               int** blkColInds, int** blkColOffsets,
                               int* blkColLens, double** remEssBCCoefs);

   //functions for getting/setting matrix or vector pointers.

   //getMatrixPtr:
   //obtain a pointer to the 'A' matrix. This should be considered a
   //constant pointer -- i.e., this class remains responsible for the
   //matrix (e.g., de-allocation upon destruction). 
   int getMatrixPtr(Data& data);

   //copyInMatrix:
   //replaces the internal matrix with a copy of the input argument, scaled
   //by the coefficient 'scalar'.

   int copyInMatrix(double scalar, const Data& data);

   //copyOutMatrix:
   //passes out a copy of the internal matrix, scaled by the coefficient
   //'scalar'.

   int copyOutMatrix(double scalar, Data& data);

   //sumInMatrix:
   //accumulate (sum) a copy of the input argument into the internal
   //matrix, scaling the input by the coefficient 'scalar'.

   int sumInMatrix(double scalar, const Data& data);

   //get/setRHSVectorPtr:
   //the same semantics apply here as for the matrixPtr functions above.

   int getRHSVectorPtr(Data& data);

   //copyInRHSVector/copyOutRHSVector/sumInRHSVector:
   //the same semantics apply here as for the matrix functions above.

   int copyInRHSVector(double scalar, const Data& data);
   int copyOutRHSVector(double scalar, Data& data);
   int sumInRHSVector(double scalar, const Data& data);

   //destroyMatrixData/destroyVectorData:
   //Utility function for destroying the matrix (or vector) in Data

   int destroyMatrixData(Data& data);
   int destroyVectorData(Data& data);

   //functions for managing multiple rhs vectors
   int setNumRHSVectors(int numRHSs, const int* rhsIDs);

   //int setRHSID:
   //set the 'current' rhs context, assuming multiple rhs vectors.
   int setRHSID(int rhsID);

   //int putInitialGuess:
   //function for setting (a subset of) the initial-guess
   //solution values (i.e., in the 'x' vector).

   int putInitialGuess(const int* eqnNumbers, const double* values,
                        int len);

   //function for getting all of the answers ('x' vector).
   int getSolution(double* answers, int len);

   //function for getting the (single) entry at equation
   //number 'eqnNumber'.
   int getSolnEntry(int eqnNumber, double& answer);

   //function for obtaining the residual vector (using current solution or
   //initial guess)
   int formResidual(double* values, int len);

   //function for launching the linear solver
   int launchSolver(int& solveStatus, int& iterations);

   int putNodalFieldData(int, int, int*, int, const double*)
     { return(0); }

   int writeSystem(const char* name);

   double* getMatrixBeginPointer() { return A_ptr_->getBeginPointer(); }

   int getMatrixOffset(int row, int col) { return A_ptr_->getOffset(row,col); }

 private:        //functions

   int createMiscStuff();

   int allocateMatrix(int** ptColIndices, int* ptRowLengths,
                       int** blkColIndices, int* blkRowLengths,
                       int* ptRowsPerBlkRow);

   int VBRmatPlusScaledMat(AztecDVBR_Matrix* A, double scalar,
                            AztecDVBR_Matrix* source);

   int MSRmatPlusScaledMat(AztecDMSR_Matrix* A, double scalar,
                            AztecDMSR_Matrix* source);

   int createBlockMatrix(int** blkColIndices,
                          int* blkRowLengths,
                          int* ptRowsPerBlkRow);

   int sumIntoBlockRow(int numBlkRows, const int* blkRows,
                        int numBlkCols, const int* blkCols,
                        const double* const* values,
                        int numPtCols,
           bool overwriteInsteadOfAccumulate);

   int copyBlockRow(int i, const int* blkRows,
                     int numBlkCols, const int* blkCols,
                     const double* const* values,
                     double* coefs);

   int modifyRHSforBCs();

   int explicitlySetDirichletBCs();

   int blockRowToPointRow(int blkRow);

   int getBlockRow(int blkRow, double*& val, int& valLen,
                   int*& blkColInds, int& blkColIndLen,
                   int& numNzBlks, int& numNNZ);

   int getBlkEqnsAndOffsets(int* ptEqns, int* blkEqns, int* blkOffsets,
                             int numEqns);

   int getBlockSize(int blkInd);

   int sumIntoPointRow(int numPtRows, const int* ptRows,
           int numPtCols, const int* ptColIndices,
           const double* const* values,
           bool overwriteInsteadOfAccumulate);

   int sumPointIntoBlockRow(int blkRow, int rowOffset,
          int blkCol, int colOffset, double value);

   int setMatrixType(const char* name);
   int selectSolver(const char* name);
   int selectPreconditioner(const char* name);
   void setSubdomainSolve(const char* name);
   void setScalingOption(const char* param);
   void setConvTest(const char* param);
   void setPreCalc(const char* param);
   void setTypeOverlap(const char* param);
   void setOverlap(const char* param);
   void setOrthog(const char* param);
   void setAuxVec(const char* param);
   void setAZ_output(const char* param);

   void recordUserParams();

   void checkForParam(const char* paramName, int numParams_,
                      char** paramStrings,
                      double& param);

   void recordUserOptions();

   void checkForOption(const char* paramName, int numParams_,
                       char** paramStrings,
                       int& param);

   int blkRowEssBCMod(int blkEqn, int blkOffset, double* val, int* blkCols,
                        int numCols, int numPtNNZ, double alpha, double gamma);

   int blkColEssBCMod(int blkRow, int blkEqn, int blkOffset, double* val,
                      int* blkCols,
                      int numCols, int numPtNNZ, double alpha, double gamma);

   void setDebugOutput(const char* path, const char* name);

   void debugOutput(const char* msg) const;

   int writeA(const char* name);
   int writeVec(Aztec_LSVector* v, const char* name);

   int messageAbort(const char* msg) const;

 private:            //variables

   MPI_Comm comm_;

   Lookup* lookup_;
   bool haveLookup_;

   int numProcs_;
   int thisProc_;
   int masterProc_;

   int* update_;
   fei::SharedPtr<Aztec_Map> map_;
   AztecDMSR_Matrix *A_;
   AztecDMSR_Matrix *A_ptr_;
   Aztec_LSVector *x_, **b_, *bc_;
   int* essBCindices_;
   int numEssBCs_;
   bool bcsLoaded_;
   bool explicitDirichletBCs_;
   bool BCenforcement_no_column_mod_;
   Aztec_LSVector *b_ptr_;
   bool matrixAllocated_;
   bool vectorsAllocated_;
   bool blkMatrixAllocated_;
   bool matrixLoaded_;
   bool rhsLoaded_;
   bool needNewPreconditioner_;

   bool tooLateToChooseBlock_;
   bool blockMatrix_;
   fei::SharedPtr<Aztec_BlockMap> blkMap_;
   AztecDVBR_Matrix *blkA_;
   AztecDVBR_Matrix *blkA_ptr_;
   int* blkUpdate_;

   AZ_MATRIX *azA_;
   AZ_PRECOND *azP_;
   bool precondCreated_;
   AZ_SCALING *azS_;
   bool scalingCreated_;

   int *aztec_options_;
   double *aztec_params_;
   double *aztec_status_;

   double* tmp_x_;
   bool tmp_x_touched_;
   double** tmp_b_;
   double* tmp_bc_;
   bool tmp_b_allocated_;

   bool ML_Vanek_;
   int numLevels_;

   int* rhsIDs_;
   int numRHSs_;

   int currentRHS_;

   int numGlobalEqns_;
   int localOffset_;
   int numLocalEqns_;

   int numGlobalEqnBlks_;
   int localBlkOffset_;
   int numLocalEqnBlks_;
   int* localBlockSizes_;

   int numNonzeroBlocks_;

   int outputLevel_;
   int numParams_;
   char** paramStrings_;

   std::string name_;
   int debugOutput_;
   int debugFileCounter_;
   char* debugPath_;
   char* debugFileName_;
   FILE* debugFile_;

   std::map<std::string,unsigned>& named_solve_counter_;
};

}//namespace fei_trilinos

#endif