FEI Version of the Day
fei_Aztec_LinSysCore.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_trilinos_macros.hpp>
00010 #include <fei_sstream.hpp>
00011 
00012 #ifdef HAVE_FEI_AZTECOO
00013 
00014 #include <stdlib.h>
00015 #include <string.h>
00016 #include <stdio.h>
00017 #include <stdexcept>
00018 
00019 #include "fei_defs.h"
00020 #include "fei_Data.hpp"
00021 #include "fei_Lookup.hpp"
00022 #include "fei_LinearSystemCore.hpp"
00023 #include "snl_fei_Utils.hpp"
00024 #include "fei_ArrayUtils.hpp"
00025 
00026 #undef fei_file
00027 #define fei_file "fei_Aztec_LinSysCore.cpp"
00028 #include "fei_ErrMacros.hpp"
00029 
00030 #include "fei_mpi.h"
00031 
00032 #ifndef FEI_SER
00033 
00034 #define AZTEC_MPI AZTEC_MPI
00035 #define AZ_MPI AZ_MPI
00036 #ifndef MPI
00037 #define MPI MPI
00038 #endif
00039 
00040 #endif
00041 
00042 #include <az_aztec.h>
00043 
00044 #include "fei_Aztec_Map.hpp"
00045 #include "fei_Aztec_BlockMap.hpp"
00046 #include "fei_Aztec_LSVector.hpp"
00047 #include "fei_AztecDMSR_Matrix.hpp"
00048 #include "fei_AztecDVBR_Matrix.hpp"
00049 
00050 #include "fei_Aztec_LinSysCore.hpp"
00051 
00052 namespace fei_trilinos {
00053 
00054 static int azlsc_solveCounter_ = 0;
00055 static int azlsc_debugFileCounter_ = 0;
00056 
00057 static std::map<std::string, unsigned> fei_aztec_named_solve_counter;
00058 
00059 //=========CONSTRUCTOR==========================================================
00060 Aztec_LinSysCore::Aztec_LinSysCore(MPI_Comm comm)
00061  : comm_(comm),
00062    lookup_(NULL),
00063    haveLookup_(false),
00064    update_(NULL),
00065    map_(),
00066    A_(NULL),
00067    A_ptr_(NULL),
00068    x_(NULL),
00069    b_(NULL),
00070    bc_(NULL),
00071    essBCindices_(NULL),
00072    numEssBCs_(0),
00073    bcsLoaded_(false),
00074    explicitDirichletBCs_(false),
00075    BCenforcement_no_column_mod_(false),
00076    b_ptr_(NULL),
00077    matrixAllocated_(false),
00078    vectorsAllocated_(false),
00079    blkMatrixAllocated_(false),
00080    matrixLoaded_(false),
00081    rhsLoaded_(false),
00082    needNewPreconditioner_(false),
00083    tooLateToChooseBlock_(false),
00084    blockMatrix_(false),
00085    blkMap_(),
00086    blkA_(NULL),
00087    blkA_ptr_(NULL),
00088    blkUpdate_(NULL),
00089    azA_(NULL),
00090    azP_(NULL),
00091    precondCreated_(false),
00092    azS_(NULL),
00093    scalingCreated_(false),
00094    aztec_options_(NULL),
00095    aztec_params_(NULL),
00096    aztec_status_(NULL),
00097    tmp_x_(NULL),
00098    tmp_x_touched_(false),
00099    tmp_b_(NULL),
00100    tmp_bc_(NULL),
00101    tmp_b_allocated_(false),
00102    ML_Vanek_(false),
00103    numLevels_(9), //arbitrary default...
00104    rhsIDs_(NULL),
00105    numRHSs_(0),
00106    currentRHS_(-1),
00107    numGlobalEqns_(0),
00108    localOffset_(0),
00109    numLocalEqns_(0),
00110    numGlobalEqnBlks_(0),
00111    localBlkOffset_(0),
00112    numLocalEqnBlks_(0),
00113    localBlockSizes_(NULL),
00114    outputLevel_(0),
00115    numParams_(0),
00116    paramStrings_(NULL),
00117    name_("dbg"),
00118    debugOutput_(0),
00119    debugFileCounter_(0),
00120    debugPath_(NULL),
00121    debugFileName_(NULL),
00122    debugFile_(NULL),
00123    named_solve_counter_(fei_aztec_named_solve_counter)
00124 {
00125    masterProc_ = 0;
00126    numProcs_ = 1;
00127    thisProc_ = 0;
00128 #ifndef FEI_SER
00129    MPI_Comm_size(comm_, &numProcs_);
00130    MPI_Comm_rank(comm_, &thisProc_);
00131 #endif
00132 
00133    aztec_options_ = new int[AZ_OPTIONS_SIZE];
00134    aztec_params_ = new double[AZ_PARAMS_SIZE];
00135 
00136    AZ_defaults(aztec_options_, aztec_params_);
00137 
00138    aztec_status_ = new double[AZ_STATUS_SIZE];
00139    for(int i=0; i<AZ_STATUS_SIZE; i++) aztec_status_[i] = 0.0;
00140 
00141    numRHSs_ = 1;
00142    rhsIDs_ = new int[numRHSs_];
00143    rhsIDs_[0] = 0;
00144 
00145    std::map<std::string,unsigned>::iterator
00146      iter = named_solve_counter_.find(name_);
00147    if (iter == named_solve_counter_.end()) {
00148      named_solve_counter_.insert(std::make_pair(name_, 0));
00149    }
00150 }
00151 
00152 //========DESTRUCTOR============================================================
00153 Aztec_LinSysCore::~Aztec_LinSysCore() {
00154 
00155   if (blkMatrixAllocated_) {
00156     delete blkA_;
00157     delete [] blkUpdate_;
00158     delete [] localBlockSizes_;
00159     blkMatrixAllocated_ = false;
00160   }
00161   if (matrixAllocated_) {
00162     delete A_;
00163     matrixAllocated_ = false;
00164   }
00165 
00166   if (precondCreated_) {
00167     AZ_precond_destroy(&azP_);
00168     precondCreated_ = false;
00169   }
00170 
00171   if (scalingCreated_) {
00172     AZ_scaling_destroy(&azS_);
00173     scalingCreated_ = false;
00174   }
00175 
00176   if (vectorsAllocated_) delete x_;
00177   delete [] tmp_x_;
00178   delete [] tmp_bc_;
00179 
00180   for(int i=0; i<numRHSs_; i++) {
00181     if (vectorsAllocated_) delete b_[i];
00182     if (tmp_b_allocated_) delete [] tmp_b_[i];
00183   }
00184   if (vectorsAllocated_) delete [] b_;
00185   if (tmp_b_allocated_) delete [] tmp_b_;
00186   tmp_b_allocated_ = false;
00187   delete [] update_;
00188 
00189   delete [] aztec_options_;
00190   delete [] aztec_params_;
00191   delete [] aztec_status_;
00192 
00193   delete [] rhsIDs_;
00194   numRHSs_ = 0;
00195 
00196   for(int j=0; j<numParams_; j++) {
00197     delete [] paramStrings_[j];
00198   }
00199   delete [] paramStrings_;
00200   numParams_ = 0;
00201 
00202   if (debugOutput_) {
00203     debugOutput_ = 0;
00204     fclose(debugFile_);
00205     delete [] debugPath_;
00206     delete [] debugFileName_;
00207   }
00208 
00209   delete [] essBCindices_;
00210   essBCindices_ = NULL;
00211   numEssBCs_ = 0;
00212   delete bc_;
00213 }
00214 
00215 //==============================================================================
00216 LinearSystemCore* Aztec_LinSysCore::clone() {
00217    return(new Aztec_LinSysCore(comm_));
00218 }
00219 
00220 //==============================================================================
00221 int Aztec_LinSysCore::parameters(int numParams, const char*const * params) {
00222 //
00223 // this function takes parameters for setting internal things like solver
00224 // and preconditioner choice, etc.
00225 //
00226    debugOutput("parameters");
00227 
00228    if (numParams == 0 || params == NULL) {
00229       debugOutput("--- no parameters");
00230       return(0);
00231    }
00232 
00233    const char* param = NULL;
00234 
00235    snl_fei::mergeStringLists(paramStrings_, numParams_,
00236                                     params, numParams);
00237 
00238    param = snl_fei::getParamValue("outputLevel",numParams,params);
00239    if (param != NULL){
00240       sscanf(param,"%d", &outputLevel_);
00241    }
00242 
00243    param = snl_fei::getParamValue("AZ_matrix_type", numParams, params);
00244    if (param != NULL){
00245       setMatrixType(param);
00246    }
00247 
00248    param = snl_fei::getParam("EXPLICIT_BC_ENFORCEMENT", numParams, params);
00249    if (param != NULL){
00250      explicitDirichletBCs_ = true;
00251    }
00252    else {
00253      explicitDirichletBCs_ = false;
00254    }
00255 
00256   param = snl_fei::getParam("BC_ENFORCEMENT_NO_COLUMN_MOD",numParams,params);
00257   if (param != NULL){
00258     BCenforcement_no_column_mod_ = true;
00259   }
00260   else {
00261     BCenforcement_no_column_mod_ = false;
00262   }
00263 
00264    param = snl_fei::getParamValue("numLevels", numParams, params);
00265    if (param != NULL){
00266       sscanf(param,"%d", &numLevels_);
00267    }
00268 
00269    bool dbgOutputParam = false;
00270    char* dbgFileName = NULL;
00271    char* dbgPath = NULL;
00272 
00273    param = snl_fei::getParamValue("debugOutput",numParams,params);
00274    if (param != NULL){
00275       dbgOutputParam = true;
00276       dbgFileName = new char[128];
00277       dbgPath = new char[strlen(param)+1];
00278       strcpy(dbgPath, param);
00279 
00280       sprintf(dbgFileName, "AZLSC.%d.%d.%d",
00281               numProcs_, thisProc_, azlsc_debugFileCounter_);
00282    }
00283 
00284    param = snl_fei::getParamValue("name",numParams, params);
00285    if(param != NULL){
00286       name_ = param;
00287 
00288       std::map<std::string,unsigned>::iterator
00289         iter = named_solve_counter_.find(name_);
00290       if (iter == named_solve_counter_.end()) {
00291         named_solve_counter_.insert(std::make_pair(name_, 0));
00292       }
00293 
00294       if (dbgOutputParam) {
00295          if (dbgFileName != NULL) delete [] dbgFileName;
00296          dbgFileName = new char[256];
00297          sprintf(dbgFileName, "AZLSC_%s.%d.%d.%d", name_.c_str(),
00298                  numProcs_, thisProc_, azlsc_debugFileCounter_);
00299       }
00300    }
00301 
00302    if (dbgOutputParam) {
00303      if (azlsc_debugFileCounter_ == azlsc_solveCounter_) {
00304        setDebugOutput(dbgPath,dbgFileName);
00305        ++azlsc_debugFileCounter_;
00306      }
00307      delete [] dbgFileName;
00308      delete [] dbgPath;
00309    }
00310 
00311    if (debugOutput_) {
00312       fprintf(debugFile_,"--- numParams %d\n",numParams);
00313       for(int i=0; i<numParams; i++){
00314          fprintf(debugFile_,"------ paramStrings[%d]: %s\n",i,
00315                  paramStrings_[i]);
00316       }
00317    }
00318 
00319    debugOutput("leaving parameters function");
00320    return(0);
00321 }
00322 
00323 //==============================================================================
00324 int Aztec_LinSysCore::setLookup(Lookup& lookup)
00325 {
00326    lookup_ = &lookup;
00327    haveLookup_ = true;
00328    return(0);
00329 }
00330 
00331 //==============================================================================
00332 int Aztec_LinSysCore::setGlobalOffsets(int len, int* nodeOffsets,
00333                                         int* eqnOffsets, int* blkEqnOffsets)
00334 {
00335   localOffset_ = eqnOffsets[thisProc_];
00336   numLocalEqns_ = eqnOffsets[thisProc_+1]-localOffset_;
00337   numGlobalEqns_ = eqnOffsets[numProcs_];
00338 
00339   localBlkOffset_ = blkEqnOffsets[thisProc_];
00340   numLocalEqnBlks_ = blkEqnOffsets[thisProc_+1]-localBlkOffset_;
00341   numGlobalEqnBlks_ = blkEqnOffsets[numProcs_];
00342 
00343   int err = createMiscStuff();
00344   if (err) return(err);
00345 
00346   if (debugOutput_) {
00347     fprintf(debugFile_, "setGlobalNodeOffsets, len: %d\n", len);
00348     for(int i=0; i<len; i++) {
00349       fprintf(debugFile_, "   nodeOffsets[%d]: %d, eqnOffsets[%d]: %d, blkEqnOffsets[%d], %d\n", i, nodeOffsets[i], i, eqnOffsets[i], i, blkEqnOffsets[i]);
00350     }
00351     fflush(debugFile_);
00352   }
00353   return(0);
00354 }
00355 
00356 //==============================================================================
00357 int Aztec_LinSysCore::setConnectivities(GlobalID elemBlock,
00358                                   int numElements,
00359                                   int numNodesPerElem,
00360                                   const GlobalID* elemIDs,
00361                                   const int* const* connNodes)
00362 {
00363    (void) elemBlock;
00364    (void) numElements;
00365    (void) numNodesPerElem;
00366    (void) elemIDs;
00367    (void) connNodes;
00368    return(0);
00369 }
00370 
00371 //==============================================================================
00372 int Aztec_LinSysCore::setMatrixType(const char* name) {
00373 
00374    if (!strcmp(name, "AZ_VBR_MATRIX")) {
00375       if (!tooLateToChooseBlock_) {
00376          blockMatrix_ = true;
00377          debugOutput("setMatrixType: chose block matrix");
00378       }
00379       else {
00380          fei::console_out() << "Aztec_LinSysCore::setMatrixType: WARNING: Too late to choose"
00381            << " the DVBR matrix; make this choice before calling "
00382            << "setMatrixStructure. DMSR will be used." << FEI_ENDL;
00383       }
00384    }
00385    else if (!strcmp(name, "AZ_MSR_MATRIX")) {
00386       //do nothing, this is the default case.
00387    }
00388    else {
00389       if (thisProc_ == 0) {
00390          FEI_COUT << "Aztec_LinSysCore: Warning: requested matrix-type <"<<name <<"> not recognized." << FEI_ENDL;
00391          FEI_COUT << "Aztec_LinSysCore: valid matrix-types are: 'AZ_MSR_MATRIX', 'AZ_VBR_MATRIX'" << FEI_ENDL;
00392       }
00393    }
00394    return(0);
00395 }
00396 
00397 //==============================================================================
00398 int Aztec_LinSysCore::setMatrixStructure(int** ptColIndices,
00399                                           int* ptRowLengths,
00400                                           int** blkColIndices,
00401                                           int* blkRowLengths,
00402                                           int* ptRowsPerBlkRow)
00403 {
00404   if (debugOutput_) {
00405     fprintf(debugFile_, "setMatrixStructure\n");
00406     for(int i=0; i<numLocalEqnBlks_; i++) {
00407       int blkEqn = localBlkOffset_+i;
00408       fprintf(debugFile_, "   blkRow %d, ptRowsPerBlkRow %d\n",
00409               blkEqn, ptRowsPerBlkRow[i]);
00410     }
00411     fflush(debugFile_);
00412   }
00413 
00414   int err = allocateMatrix(ptColIndices, ptRowLengths, blkColIndices,
00415                            blkRowLengths, ptRowsPerBlkRow);
00416   return(err);
00417 }
00418 
00419 //==============================================================================
00420 int Aztec_LinSysCore::createMiscStuff()
00421 {
00422 //
00423 //This function is where we establish the structures/objects associated
00424 //with the linear algebra library. i.e., do initial allocations, etc.
00425 //
00426    if (debugOutput_) {
00427       fprintf(debugFile_,
00428               "createMiscStuff: numRHSs_: %d\n", numRHSs_);
00429       fflush(debugFile_);
00430    }
00431 
00432    if (numLocalEqns_ > numGlobalEqns_)
00433       messageAbort("createMiscStuff: numLocalEqns > numGlobalEqns");
00434 
00435    if (0 > localOffset_)
00436       messageAbort("createMiscStuff: localOffset_ out of range");
00437 
00438    update_ = numLocalEqns_ > 0 ? new int[numLocalEqns_] : NULL;
00439    if (numLocalEqns_ > 0 && update_ == NULL) {
00440      return(-1);
00441    }
00442 
00443    int i, j;
00444    for(j=0; j<numLocalEqns_; j++) update_[j] = localOffset_+j;
00445 
00446    azS_ = AZ_scaling_create();
00447    if (azS_ != NULL) scalingCreated_ = true;
00448    else {
00449      fei::console_out() << "Aztec_LinSysCore::createMiscStuff ERROR: failed to create scaling"
00450           << FEI_ENDL;
00451      return(-1);
00452    }
00453 
00454    if (numRHSs_ <= 0)
00455       messageAbort("numRHSs_==0. Out of scope or destructor already called?");
00456 
00457     tmp_x_ = numLocalEqns_ > 0 ? new double[numLocalEqns_] : NULL;
00458     if (numLocalEqns_ > 0 && tmp_x_ == NULL) return(-1);
00459 
00460     tmp_bc_ = numLocalEqns_ > 0 ? new double[numLocalEqns_] : NULL;
00461     if (numLocalEqns_ > 0 && tmp_bc_ == NULL) return(-1);
00462 
00463     for(i=0; i<numLocalEqns_; i++){
00464         tmp_x_[i] = 0.0;
00465         tmp_bc_[i] = 0.0;
00466     }
00467 
00468     if (!tmp_b_allocated_) {
00469        tmp_b_ = new double*[numRHSs_];
00470 
00471        for(j=0; j<numRHSs_; j++) {
00472            tmp_b_[j] = new double[numLocalEqns_];
00473            for(i=0; i<numLocalEqns_; i++) {
00474                tmp_b_[j][i] = 0.0;
00475            }
00476        }
00477 
00478        tmp_b_allocated_ = true;
00479     }
00480 
00481    if (currentRHS_ < 0) currentRHS_ = 0;
00482    return(0);
00483 }
00484 
00485 //==============================================================================
00486 int Aztec_LinSysCore::allocateMatrix(int** ptColIndices,
00487                                       int* ptRowLengths,
00488                                       int** blkColIndices,
00489                                       int* blkRowLengths,
00490                                       int* ptRowsPerBlkRow)
00491 {
00492   int err;
00493   if (blockMatrix_) {
00494     err = createBlockMatrix(blkColIndices, blkRowLengths, ptRowsPerBlkRow);
00495     return(err);
00496   }
00497 
00498   tooLateToChooseBlock_ = true;
00499 
00500   map_.reset( new Aztec_Map(numGlobalEqns_, numLocalEqns_, update_, localOffset_, comm_));
00501 
00502   A_ = new AztecDMSR_Matrix(map_);
00503 
00504   if (A_ptr_ == NULL) A_ptr_ = A_;
00505 
00506   //
00507   //we're going to have to look through the colIndices and see if any
00508   //of the rows do NOT have an entry on the diagonal. For each row that
00509   //does have an entry on the diagonal, we subtract 1 from the row-length
00510   //since the AztecDMSR_Matrix::allocate function wants the length of each
00511   //row NOT including any entry on the diagonal.
00512   //
00513 
00514   int* row_lengths = numLocalEqns_ > 0 ? new int[numLocalEqns_] : NULL;
00515 
00516   for (int i = 0; i < numLocalEqns_; i++) {
00517     if (fei::searchList(localOffset_+i,
00518                             ptColIndices[i], ptRowLengths[i]) >= 0) {
00519       row_lengths[i] = ptRowLengths[i] - 1;
00520     }
00521     else {
00522       row_lengths[i] = ptRowLengths[i];
00523     }
00524   }
00525 
00526   // so now we know all the row lengths, and can configure our matrix.
00527 
00528   A_ptr_->allocate(row_lengths, ptColIndices);
00529 
00530   delete [] row_lengths;
00531 
00532   matrixAllocated_ = true;
00533   return(0);
00534 }
00535 
00536 //==============================================================================
00537 int Aztec_LinSysCore::createBlockMatrix(int** blkColIndices,
00538                                          int* blkRowLengths,
00539                                          int* ptRowsPerBlkRow)
00540 {
00541   int i;
00542 
00543   blkUpdate_ = new int[numLocalEqnBlks_];
00544 
00545   for(i=0; i<numLocalEqnBlks_; i++) {
00546     blkUpdate_[i] = localBlkOffset_ + i;
00547   }
00548 
00549   blkMap_.reset(new Aztec_BlockMap(numGlobalEqns_, numLocalEqns_,
00550                                update_, localOffset_, comm_,
00551                                numGlobalEqnBlks_, numLocalEqnBlks_,
00552                                blkUpdate_,
00553                                localBlkOffset_, ptRowsPerBlkRow));
00554 
00555   blkA_ = new AztecDVBR_Matrix(blkMap_);
00556 
00557   if (blkA_ptr_ == NULL) blkA_ptr_ = blkA_;
00558 
00559   localBlockSizes_ = new int[numLocalEqnBlks_];
00560 
00561   //now we need to count how many non-zero blocks there are.
00562   numNonzeroBlocks_ = 0;
00563   for(i=0; i<numLocalEqnBlks_; i++) {
00564 
00565     numNonzeroBlocks_ += blkRowLengths[i];
00566 
00567     localBlockSizes_[i] = ptRowsPerBlkRow[i];
00568   }
00569 
00570   //and now we need to flatten the list of block-column-indices into a 1-D
00571   //list to give to the AztecDVBR_Matrix::allocate function.
00572   int* blk_col_inds = new int[numNonzeroBlocks_];
00573 
00574   int offset = 0;
00575   for(i=0; i<numLocalEqnBlks_; i++) {
00576     for(int j=0; j<blkRowLengths[i]; j++) {
00577       blk_col_inds[offset++] = blkColIndices[i][j];
00578     }
00579   }
00580 
00581   //finally we're ready to allocate our matrix.
00582   blkA_->allocate(blkRowLengths, blk_col_inds);
00583 
00584   delete [] blk_col_inds;
00585 
00586   blkMatrixAllocated_ = true;
00587   return(0);
00588 }
00589 
00590 //==============================================================================
00591 int Aztec_LinSysCore::resetMatrixAndVector(double s)
00592 {
00593    if (blkMatrixAllocated_) blkA_ptr_->put(s);
00594    if (matrixAllocated_) A_ptr_->put(s);
00595 
00596    if (rhsLoaded_) {
00597       for(int i=0; i<numRHSs_; i++){
00598          b_[i]->put(s);
00599       }
00600    }
00601 
00602    if (b_ptr_ != NULL) b_ptr_->put(s);
00603 
00604    if (bc_ != NULL) bc_->put(s);
00605 
00606    if (!tmp_b_allocated_) {
00607       for(int j=0; j<numRHSs_; j++) {
00608          tmp_b_[j] = new double[numLocalEqns_];
00609       }
00610    }
00611    matrixLoaded_ = false;
00612    rhsLoaded_ = false;
00613    bcsLoaded_ = false;
00614 
00615    delete [] tmp_bc_;
00616    tmp_bc_ = new double[numLocalEqns_];
00617 
00618    for(int ii=0; ii<numLocalEqns_; ii++) tmp_bc_[ii] = s;
00619 
00620    delete [] essBCindices_;
00621    essBCindices_ = NULL;
00622    numEssBCs_ = 0;
00623 
00624    for(int j=0; j<numRHSs_; j++) {
00625       for(int i=0; i<numLocalEqns_; i++) {
00626          tmp_b_[j][i] = s;
00627       }
00628    }
00629    return(0);
00630 }
00631 
00632 //==============================================================================
00633 int Aztec_LinSysCore::resetMatrix(double s) {
00634    if (blkMatrixAllocated_) blkA_ptr_->put(s);
00635    if (matrixAllocated_) A_ptr_->put(s);
00636 
00637    matrixLoaded_ = false;
00638    bcsLoaded_ = false;
00639 
00640    if (bc_ != NULL) bc_->put(s);
00641 
00642    delete [] tmp_bc_;
00643    tmp_bc_ = new double[numLocalEqns_];
00644 
00645    for(int ii=0; ii<numLocalEqns_; ii++) tmp_bc_[ii] = s;
00646 
00647    delete [] essBCindices_;
00648    essBCindices_ = NULL;
00649    numEssBCs_ = 0;
00650 
00651    return(0);
00652 }
00653 
00654 //==============================================================================
00655 int Aztec_LinSysCore::resetRHSVector(double s) {
00656 
00657    if (rhsLoaded_) {
00658       for(int i=0; i<numRHSs_; i++){
00659          b_[i]->put(s);
00660       }
00661    }
00662 
00663    if (b_ptr_ != NULL) b_ptr_->put(s);
00664 
00665    if (!tmp_b_allocated_) {
00666       for(int j=0; j<numRHSs_; j++) {
00667          tmp_b_[j] = new double[numLocalEqns_];
00668       }
00669    }
00670    rhsLoaded_ = false;
00671    bcsLoaded_ = false;
00672 
00673    for(int j=0; j<numRHSs_; j++) {
00674      double* cur_b = tmp_b_[j];
00675      for(int i=0; i<numLocalEqns_; i++) {
00676        cur_b[i] = s;
00677      }
00678    }
00679    return(0);
00680 }
00681 
00682 //==============================================================================
00683 int Aztec_LinSysCore::sumIntoSystemMatrix(int numPtRows, const int* ptRows,
00684                             int numPtCols, const int* ptCols,
00685                             const double* const* values)
00686 {
00687   matrixLoaded_ = false;
00688 
00689   return(sumIntoPointRow(numPtRows, ptRows, numPtCols, ptCols, values, false));
00690 }
00691 
00692 //==============================================================================
00693 int Aztec_LinSysCore::sumIntoSystemMatrix(int numPtRows, const int* ptRows,
00694                                            int numPtCols, const int* ptCols,
00695                                            int numBlkRows, const int* blkRows,
00696                                            int numBlkCols, const int* blkCols,
00697                                            const double* const* values)
00698 {
00699    matrixLoaded_ = false;
00700 
00701    if ((A_ptr_ == NULL) && (blkA_ptr_ == NULL))
00702       messageAbort("sumIntoSystemMatrix: matrix is NULL.");
00703 
00704    if (blockMatrix_) {
00705       return(sumIntoBlockRow(numBlkRows, blkRows, numBlkCols, blkCols,
00706                       values, numPtCols, false));
00707    }
00708    else {
00709       return( sumIntoPointRow(numPtRows, ptRows, numPtCols, ptCols, values,
00710                               false));
00711    }
00712 }
00713 
00714 //==============================================================================
00715 int Aztec_LinSysCore::putIntoSystemMatrix(int numPtRows, const int* ptRows,
00716                                            int numPtCols, const int* ptCols,
00717                                            const double* const* values)
00718 {
00719   matrixLoaded_ = false;
00720 
00721   return(sumIntoPointRow(numPtRows, ptRows, numPtCols, ptCols, values, true));
00722 }
00723 
00724 //==============================================================================
00725 int Aztec_LinSysCore::getMatrixRowLength(int row, int& length)
00726 {
00727   if (!blockMatrix_) {
00728     if (row >= localOffset_ && row < (localOffset_+numLocalEqns_)) {
00729       length = A_ptr_->rowLength(row);
00730       return(0);
00731     }
00732   }
00733   else {
00734     if (!haveLookup_) return(-1);
00735     int blkRow = lookup_->ptEqnToBlkEqn(row);
00736     if (blkRow < 0) return(-1);
00737     return( blkA_ptr_->getNumNonzerosPerRow(blkRow, length) );
00738   }
00739 
00740   return(-1);
00741 }
00742 
00743 //==============================================================================
00744 int Aztec_LinSysCore::getMatrixRow(int row, double* coefs,
00745                                    int* indices,
00746                                    int len, int& rowLength)
00747 {
00748   int err = getMatrixRowLength(row, rowLength);
00749   if (err != 0) return(-1);
00750 
00751   int* tmpIndices = indices;
00752   double* tmpCoefs = coefs;
00753 
00754   if (len < rowLength) {
00755     tmpIndices = new int[rowLength];
00756     tmpCoefs = new double[rowLength];
00757   }
00758 
00759   if (!blockMatrix_) {
00760     A_ptr_->getRow(row, rowLength, tmpCoefs, tmpIndices);
00761   }
00762   else {
00763     fei::console_out() << "Aztec_LinSysCore::getMatrixRow: not implemented." << FEI_ENDL;
00764     return(-1);
00765   }
00766 
00767   if (len < rowLength) {
00768     for(int i=0; i<len; i++) {
00769       indices[i] = tmpIndices[i];
00770       coefs[i] = tmpCoefs[i];
00771     }
00772     delete [] tmpIndices;
00773     delete [] tmpCoefs;
00774   }
00775 
00776   return(0);
00777 }
00778 
00779 //==============================================================================
00780 int Aztec_LinSysCore::sumIntoBlockRow(int numBlkRows, const int* blkRows,
00781                                        int numBlkCols, const int* blkCols,
00782                                        const double* const* values,
00783                                        int numPtCols,
00784                                       bool overwriteInsteadOfAccumulate)
00785 {
00786   int i;
00787   //first, let's figure out which of the incoming blkRows is the biggest --
00788   //i.e., which contains the most point-rows.
00789   int maxBlkSize = 0;
00790 
00791   for(i=0; i<numBlkRows; i++) {
00792     int thisSize = lookup_->getBlkEqnSize(blkRows[i]);
00793 
00794     if (maxBlkSize < thisSize) maxBlkSize = thisSize;
00795   }
00796 
00797   //now we can allocate an array to hold the values for passing each block
00798   //row into the block-matrix.
00799   double* coefs = new double[numPtCols*maxBlkSize];
00800 
00801   int rowOffset = 0;
00802   for(i=0; i<numBlkRows; i++) {
00803     //now, copy the values for this block row into the coefs array.
00804     copyBlockRow(i, blkRows, numBlkCols, blkCols, &(values[rowOffset]), coefs);
00805 
00806     //now shove it into the matrix.
00807     if (overwriteInsteadOfAccumulate) {
00808       int err = blkA_ptr_->putBlockRow(blkRows[i], coefs,
00809                                        (int*)blkCols, numBlkCols);
00810       if (err != 0) {
00811         fei::console_out() << thisProc_ << " DVBR::putBlockRow failed." << FEI_ENDL;
00812         return(err);
00813       }
00814     }
00815     else {
00816       int err = blkA_ptr_->sumIntoBlockRow(blkRows[i], coefs,
00817                                            (int*)blkCols, numBlkCols);
00818       if (err != 0) {
00819         fei::console_out() << thisProc_ << " DVBR::sumIntoBlockRow failed." << FEI_ENDL;
00820         return(err);
00821       }
00822     }
00823 
00824     rowOffset += lookup_->getBlkEqnSize(blkRows[i]);
00825   }
00826 
00827   delete [] coefs;
00828   return(0);
00829 }
00830 
00831 //==============================================================================
00832 int Aztec_LinSysCore::copyBlockRow(int i, const int* blkRows,
00833                                     int numBlkCols, const int* blkCols,
00834                                     const double* const* values,
00835                                     double* coefs){
00836 
00837    int rowSize = 0, colSize = 0;
00838    //coefs needs to contain the entries for each block together, and those
00839    //entries need to be in column-major order
00840    int colOffset = 0;
00841    int coefOffset = 0;
00842    for(int b=0; b<numBlkCols; b++) {
00843 
00844       int err = blkA_ptr_->getBlockSize(blkRows[i], blkCols[b],
00845                                             rowSize, colSize);
00846       if (err) {
00847          fei::console_out() << "Aztec_LSC:copyBlockRow: ERROR in getBlockSize" << FEI_ENDL;
00848          return(-1);
00849       }
00850 
00851       for(int j=colOffset; j<colOffset+colSize; j++) {
00852          for(int r=0; r<rowSize; r++) {
00853             coefs[coefOffset++] = values[r][j];
00854          }
00855       }
00856       colOffset += colSize;
00857    }
00858    return(0);
00859 }
00860 
00861 //==============================================================================
00862 int Aztec_LinSysCore::getBlockSize(int blkInd) {
00863 
00864    int localBlkRow = blkInd - localBlkOffset_;
00865    if (localBlkRow >= 0 && localBlkRow < numLocalEqnBlks_) {
00866       //if this is a local blkInd, we have its size in an array.
00867       return(localBlockSizes_[localBlkRow]);
00868    }
00869    else {
00870       //else it ain't local, so we'll get its size from the matrix because
00871       //we know that the matrix obtained it when it allocated itself.
00872 
00873       int numRemoteBlocks = blkA_ptr_->getNumRemoteBlocks();
00874       int* remoteInds = blkA_ptr_->getRemoteBlockIndices();
00875       int* remoteSizes = blkA_ptr_->getRemoteBlockSizes();
00876 
00877       int ins;
00878       int index = fei::binarySearch(blkInd, remoteInds, numRemoteBlocks, ins);
00879       if (index < 0)
00880          messageAbort("getBlockSize: can't find blkInd.");
00881 
00882       return(remoteSizes[index]);
00883    }
00884 }
00885 
00886 //==============================================================================
00887 int Aztec_LinSysCore::sumIntoPointRow(int numPtRows, const int* ptRows,
00888                                        int numPtCols, const int* ptColIndices,
00889                                        const double* const* values,
00890                                       bool overwriteInsteadOfAccumulate)
00891 {
00892   int i, j;
00893   if (debugOutput_) {
00894     fprintf(debugFile_,"sumIntoPointRow, %d rows\n", numPtRows);
00895     for(i=0; i<numPtRows; ++i) {
00896       for(j=0; j<numPtCols; ++j) {
00897         fprintf(debugFile_,"  sipr row %d, col %d, value: %e\n", ptRows[i],
00898                 ptColIndices[j], values[i][j]);
00899       }
00900     }
00901     fflush(debugFile_);
00902   }
00903 
00904   if (!blockMatrix_) {
00905     if (overwriteInsteadOfAccumulate) {
00906       for(i=0; i<numPtRows; ++i) {
00907         CHK_ERR( A_ptr_->putRow(ptRows[i], numPtCols, values[i], ptColIndices) );
00908       }
00909     }
00910     else {
00911       int err = A_ptr_->sumIntoRow(numPtRows, ptRows, numPtCols, ptColIndices,
00912                                    values);
00913       if (err != 0) {
00914         FEI_OSTRINGSTREAM osstr;
00915         osstr << "Aztec_LinSysCore::sumIntoPointRow ERROR calling A_ptr->sumIntoRow";
00916         throw std::runtime_error(osstr.str());
00917       }
00918     }
00919 
00920     return(0);
00921   }
00922 
00923   if (!haveLookup_) {
00924     messageAbort("sumIntoPointRow: need lookup object, don't have it.");
00925   }
00926 
00927   int* blkIntData = new int[numPtRows*2 + numPtCols*2];
00928   int* blkRows = blkIntData;
00929   int* blkRowOffsets = blkIntData+numPtRows;
00930   int* blkCols = blkIntData+2*numPtRows;
00931   int* blkColOffsets = blkIntData+2*numPtRows+numPtCols;;
00932 
00933   //now fill the blkRows and blkCols lists.
00934 
00935   for(i=0; i<numPtRows; i++) {
00936     blkRows[i] = lookup_->ptEqnToBlkEqn(ptRows[i]);
00937     blkRowOffsets[i] = lookup_->getOffsetIntoBlkEqn(blkRows[i], ptRows[i]);
00938   }
00939 
00940   for(i=0; i<numPtCols; i++) {
00941     blkCols[i] = lookup_->ptEqnToBlkEqn(ptColIndices[i]);
00942     if (blkCols[i] < 0) {
00943       fei::console_out() << thisProc_ << " lookup ptEqnToBlkEqn("<<ptColIndices[i]<<"): "
00944            << blkCols[i] << FEI_ENDL;
00945       messageAbort("negative blk-col");
00946     }
00947 
00948     blkColOffsets[i] = lookup_->getOffsetIntoBlkEqn(blkCols[i], ptColIndices[i]);
00949   }
00950 
00951   for(i=0; i<numPtRows; i++) {
00952     for(j=0; j<numPtCols; j++) {
00953       sumPointIntoBlockRow(blkRows[i], blkRowOffsets[i],
00954                            blkCols[j], blkColOffsets[j], values[i][j]);
00955     }
00956   }
00957 
00958   delete [] blkIntData;
00959   return(0);
00960 }
00961 
00962 //==============================================================================
00963 int Aztec_LinSysCore::sumPointIntoBlockRow(int blkRow, int rowOffset,
00964                                             int blkCol, int colOffset,
00965                                             double value)
00966 {
00967    int rowSize = lookup_->getBlkEqnSize(blkRow);
00968    int colSize = lookup_->getBlkEqnSize(blkCol);
00969 
00970    int len = rowSize*colSize;
00971    if (len <= 0) {
00972       fei::console_out() << thisProc_ << ", ALSC::spibr: blkRow: " << blkRow << ", blkCol: " << blkCol << ", rowSize: " << rowSize << ", colSize: " << colSize << FEI_ENDL;
00973       messageAbort("sumPointIntoBlockRow: len <= 0");
00974    }
00975 
00976    double* val = new double[rowSize*colSize];
00977 
00978    for(int i=0; i<len; i++) val[i] = 0.0;
00979 
00980    val[colOffset*rowSize + rowOffset] = value;
00981 
00982    int err = blkA_ptr_->sumIntoBlockRow(blkRow, val, &blkCol, 1);
00983    if (err != 0) {
00984      fei::console_out() << thisProc_ << " DVBR::sumIntoBlockRow failed" << FEI_ENDL;
00985      return(err);
00986    }
00987 
00988    delete [] val;
00989    return(0);
00990 }
00991 
00992 //==============================================================================
00993 int Aztec_LinSysCore::sumIntoRHSVector(int num,
00994                                        const double* values,
00995                                        const int* indices)
00996 {
00997   //
00998   //This function scatters (accumulates) values into the linear-system's
00999   //currently selected RHS vector.
01000   //
01001   // num is how many values are being passed,
01002   // indices holds the global 'row-numbers' into which the values go,
01003   // and values holds the actual coefficients to be scattered.
01004   //
01005   rhsLoaded_ = false;
01006   double* cur_b = tmp_b_[currentRHS_];
01007 
01008   if (debugOutput_) {
01009     for(int i=0; i<num; ++i) {
01010       fprintf(debugFile_,"sumIntoRHS %d, %e\n", indices[i], values[i]);
01011       fflush(debugFile_);
01012     }
01013   }
01014 
01015   for(int i=0; i<num; i++){
01016     int localRow = indices[i] - localOffset_;
01017     if (localRow < 0 || localRow > numLocalEqns_) continue;
01018 
01019     cur_b[localRow] += values[i];
01020   }
01021   return(0);
01022 }
01023 
01024 //==============================================================================
01025 int Aztec_LinSysCore::putIntoRHSVector(int num, const double* values,
01026                                        const int* indices)
01027 {
01028 //
01029 //This function scatters (puts) values into the linear-system's
01030 //currently selected RHS vector.
01031 //
01032 // num is how many values are being passed,
01033 // indices holds the global 'row-numbers' into which the values go,
01034 // and values holds the actual coefficients to be scattered.
01035 //
01036    rhsLoaded_ = false;
01037 
01038    for(int i=0; i<num; i++){
01039      int localRow = indices[i] - localOffset_;
01040      if (localRow < 0 || localRow > numLocalEqns_) continue;
01041 
01042      if (debugOutput_) {
01043        fprintf(debugFile_,"putIntoRHS %d, %e\n", indices[i], values[i]);
01044        fflush(debugFile_);
01045      }
01046 
01047      tmp_b_[currentRHS_][localRow] = values[i];
01048    }
01049    return(0);
01050 }
01051 
01052 //==============================================================================
01053 int Aztec_LinSysCore::getFromRHSVector(int num, double* values,
01054                                        const int* indices)
01055 {
01056   //
01057   //This function retrieves values from the linear-system's
01058   //currently selected RHS vector.
01059   //
01060   // num is how many values are being requested,
01061   // indices holds the global 'row-numbers' for which values are requested,
01062   // and values holds the actual coefficients to be returned.
01063   //
01064   //if non-local indices are supplied, the corresponding positions in the values
01065   //array are not referenced.
01066   //
01067 
01068   for(int i=0; i<num; i++){
01069     int localRow = indices[i] - localOffset_;
01070     if (localRow < 0 || localRow > numLocalEqns_) continue;
01071 
01072     values[i] = tmp_b_[currentRHS_][localRow];
01073   }
01074   return(0);
01075 }
01076 
01077 //==============================================================================
01078 int Aztec_LinSysCore::matrixLoadComplete() {
01079 
01080    if (debugOutput_) {
01081       fprintf(debugFile_,"matrixLoadComplete\n");
01082       fflush(debugFile_);
01083    }
01084 
01085    if (matrixLoaded_ && rhsLoaded_) return(0);
01086 
01087    fei::SharedPtr<Aztec_Map> tmpMap;
01088    int* data_org = NULL;
01089 
01090    if (blockMatrix_) {
01091       tmpMap = blkMap_;
01092    }
01093    else {
01094       tmpMap = map_;
01095    }
01096 
01097    if (!matrixLoaded_) {
01098      if (blockMatrix_) {
01099        if (!blkA_ptr_->isLoaded()) blkA_ptr_->loadComplete();
01100      }
01101      else {
01102        if (!A_ptr_->isFilled()) {
01103          A_ptr_->fillComplete();
01104        }
01105      }
01106 
01107      matrixLoaded_ = true;
01108 
01109      needNewPreconditioner_ = true;
01110    }
01111 
01112    if (blockMatrix_) data_org = blkA_ptr_->getData_org();
01113    else data_org = A_ptr_->getAZ_MATRIX_PTR()->data_org;
01114 
01115    Aztec_LSVector* tmp = NULL;
01116 
01117    //if x_ is not null, then it probably has a previous solution in it that
01118    //we might as well keep for the next initial guess unless the user has
01119    //specifically set an initial guess.
01120    if (x_ != NULL) {
01121       tmp = new Aztec_LSVector(*x_);
01122       *tmp = *x_;
01123    }
01124 
01125    //if x_ hasn't been allocated yet, we better do that now.
01126    if (x_ == NULL) x_ = new Aztec_LSVector(tmpMap, data_org);
01127    if (bc_ == NULL) bc_ = new Aztec_LSVector(tmpMap, data_org);
01128 
01129    //if we did save off a copy of x_ above, let's put it back now.
01130    if (tmp != NULL) {
01131       *x_ = *tmp;
01132       delete tmp;
01133    }
01134 
01135    //if we've put boundary-condition data in tmp_bc_ then we better move it
01136    //into bc_ now.
01137    if (tmp_bc_ != NULL) {
01138      for(int j=0; j<numEssBCs_; j++) {
01139        int index = essBCindices_[j];
01140        (*bc_)[index] = tmp_bc_[index-localOffset_];
01141      }
01142    }
01143 
01144    if (tmp_x_touched_) {
01145       //and now we can fill x_ with the stuff we've been holding in
01146       //tmp_x_, if tmp_x_ has been touched... i.e., if the user loaded an
01147       //initial guess into it.
01148       for(int i=0; i<numLocalEqns_; i++){
01149          (*x_)[i+localOffset_] = tmp_x_[i];
01150       }
01151    }
01152 
01153    //now we're going to get the AZ_MATRIX ptr out of the AztecDMSR matrix
01154    //wrapper class.
01155 
01156    if (blockMatrix_) azA_ = blkA_ptr_->getAZ_MATRIX_Ptr();
01157    else azA_ = A_ptr_->getAZ_MATRIX_PTR();
01158 
01159    if (rhsLoaded_) return(0);
01160 
01161    if (b_ != NULL) {
01162       for(int i=0; i<numRHSs_; i++) delete b_[i];
01163       delete [] b_;
01164    }
01165 
01166    b_ = new Aztec_LSVector*[numRHSs_];
01167    for(int j=0; j<numRHSs_; j++) {
01168       b_[j] = new Aztec_LSVector(tmpMap, data_org);
01169       if (b_[j] == NULL) return(-1);
01170 
01171       //now fill b_[j] with the stuff we've been holding in tmp_b_[j].
01172       for(int i=0; i<numLocalEqns_; i++){
01173          (*(b_[j]))[i+localOffset_] = tmp_b_[j][i];
01174       }
01175    }
01176 
01177    b_ptr_ = b_[currentRHS_];
01178    vectorsAllocated_ = true;
01179 
01180    if (!bcsLoaded_) modifyRHSforBCs();
01181    bcsLoaded_ = false;
01182 
01183    rhsLoaded_ = true;
01184 
01185    if (debugOutput_) {
01186       fprintf(debugFile_,"leaving matrixLoadComplete\n");
01187       fflush(debugFile_);
01188    }
01189    return(0);
01190 }
01191 
01192 //==============================================================================
01193 int Aztec_LinSysCore::getBlkEqnsAndOffsets(int* ptEqns, int* blkEqns,
01194                                             int* blkOffsets, int numEqns)
01195 {
01196    if (!haveLookup_) messageAbort("getBlkEqnsAndOffsets: don't have Lookup");
01197 
01198    for(int i=0; i<numEqns; i++) {
01199       blkEqns[i] = lookup_->ptEqnToBlkEqn(ptEqns[i]);
01200       if (blkEqns[i] < 0) {
01201         fei::console_out() << thisProc_ << "ptEqn: " << ptEqns[i] << ", blkEqn: " << blkEqns[i]
01202              << FEI_ENDL;
01203         messageAbort("getBlkEqnsAndOffsets: blk-eqn lookup failed.");
01204       }
01205 
01206       blkOffsets[i] = lookup_->getOffsetIntoBlkEqn(blkEqns[i], ptEqns[i]);
01207       if (blkOffsets[i] < 0) {
01208          messageAbort("getBlkEqnsAndOffsets: blk-offset lookup failed.");
01209       }
01210    }
01211    return(0);
01212 }
01213 
01214 //==============================================================================
01215 int Aztec_LinSysCore::enforceEssentialBC(int* globalEqn,
01216                                          double* alpha,
01217                                          double* gamma, int len)
01218 {
01219   //
01220   //This function must enforce an essential boundary condition on each local
01221   //equation in 'globalEqn'. This means that the following modification
01222   //should be made to A and b, for each globalEqn:
01223   //
01224   //for(each locally-owned equation i in the globalEqn array){
01225   //   zero row i and put 1.0 on the diagonal
01226   //   for(each row j in column i) {
01227   //      if (i==j) b[i] = gamma/alpha;
01228   //      else b[j] -= (gamma/alpha)*A[j,i];
01229   //      A[j,i] = 0.0;
01230   //   }
01231   //}
01232   //
01233 
01234   if (len == 0) return(0);
01235 
01236   std::vector<int> bcEqns; bcEqns.reserve(len);
01237   std::vector<int> indirect; indirect.reserve(len);
01238   for(int i=0; i<len; ++i) {
01239     bcEqns.push_back(globalEqn[i]);
01240     indirect.push_back(i);
01241   }
01242 
01243   fei::insertion_sort_with_companions<int>(len, &bcEqns[0], &indirect[0]);
01244 
01245   if (debugOutput_) {
01246     fprintf(debugFile_,"numEssBCs: %d\n", len);
01247     for(int ii=0; ii<len; ++ii) {
01248       fprintf(debugFile_, "   EssBC eqn %d, alpha %e gamma %e\n",
01249               bcEqns[ii], alpha[indirect[ii]], gamma[indirect[ii]]);
01250     }
01251     fflush(debugFile_);
01252   }
01253 
01254   bcsLoaded_ = true;
01255 
01256   int localEnd = localOffset_ + numLocalEqns_ - 1;
01257 
01258   if (debugOutput_) {
01259     fprintf(debugFile_,"localOffset_: %d, localEnd: %d\n", localOffset_, localEnd);
01260     fflush(debugFile_);
01261   }
01262 
01263   {
01264     int* newBCindices = new int[numEssBCs_+len];
01265     int ii;
01266     for(ii=0; ii<numEssBCs_; ii++) newBCindices[ii] = essBCindices_[ii];
01267     int offset = 0;
01268     for(ii=0; ii<len; ii++) {
01269       if ((localOffset_ <= globalEqn[ii]) && (globalEqn[ii] <= localEnd)){
01270         newBCindices[numEssBCs_+offset++] = globalEqn[ii];
01271       }
01272     }
01273 
01274     if (offset > 0) {
01275       delete [] essBCindices_;
01276       essBCindices_ = newBCindices;
01277       numEssBCs_ += offset;
01278     }
01279   }
01280 
01281   if (blockMatrix_) {
01282     int* blkIntData = new int[len*2];
01283     int* blkEqns = blkIntData;
01284     int* blkOffsets = blkIntData+len;
01285 
01286     getBlkEqnsAndOffsets(globalEqn, blkEqns, blkOffsets, len);
01287 
01288     enforceBlkEssentialBC(blkEqns, blkOffsets, alpha, gamma, len);
01289 
01290     delete [] blkIntData;
01291 
01292     return(0);
01293   }
01294 
01295   for(int i=0; i<len; i++) {
01296 
01297     int globalEqn_i = bcEqns[i];
01298 
01299     //if globalEqn[i] isn't local, then the processor that owns it
01300     //should be running this code too. Otherwise there's trouble...
01301 
01302     if ((localOffset_ > globalEqn_i) || (globalEqn_i > localEnd)) continue;
01303 
01304     //zero this row, except for the diagonal coefficient.
01305     A_ptr_->setDiagEntry(globalEqn_i, 1.0);
01306     int* offDiagIndices = NULL;
01307     double* offDiagCoefs = NULL;
01308     int offDiagLength = 0;
01309     A_ptr_->getOffDiagRowPointers(globalEqn_i, offDiagIndices,
01310                                   offDiagCoefs, offDiagLength);
01311 
01312     for(int jjj=0; jjj<offDiagLength; jjj++) offDiagCoefs[jjj] = 0.0;
01313 
01314     //also, make the rhs modification here.
01315     double bc_term = gamma[indirect[i]]/alpha[indirect[i]];
01316     double rhs_term = bc_term;
01317     if (explicitDirichletBCs_) rhs_term = 0.0;
01318 
01319     if (rhsLoaded_) {
01320       (*b_ptr_)[globalEqn_i] = rhs_term;
01321       (*bc_)[globalEqn_i] = bc_term;
01322     }
01323     else {
01324       tmp_b_[currentRHS_][globalEqn_i -localOffset_] = rhs_term;
01325       tmp_bc_[globalEqn_i -localOffset_] = bc_term;
01326     }
01327   }
01328 
01329   if (BCenforcement_no_column_mod_ == true) {
01330     return(0);
01331   }
01332 
01333   for(int row=localOffset_; row<=localEnd; row++) {
01334 
01335     int insertPoint = -1;
01336     int index = fei::binarySearch(row, &bcEqns[0], len, insertPoint);
01337     if (index >= 0) continue;
01338 
01339     int* offDiagIndices2 = NULL;
01340     double* offDiagCoefs2 = NULL;
01341     int offDiagLength2 = 0;
01342     A_ptr_->getOffDiagRowPointers(row, offDiagIndices2,
01343                                   offDiagCoefs2, offDiagLength2);
01344 
01345     //look through this row to find the non-zeros in positions that
01346     //correspond to eqns in bcEqns and make the appropriate modification.
01347 
01348     for(int j=0; j<offDiagLength2; j++) {
01349 
01350       int col_index = A_ptr_->getAztec_Map()->getTransformedEqn(offDiagIndices2[j]);
01351 
01352       int idx = fei::binarySearch(col_index, &bcEqns[0], len, insertPoint);
01353       if (idx < 0) continue;
01354 
01355       double bc_term = gamma[indirect[idx]]/alpha[indirect[idx]];
01356 
01357       double value = offDiagCoefs2[j]*bc_term;
01358 
01359       if (rhsLoaded_) {
01360         (*b_ptr_)[row] -= value;
01361         (*bc_)[row] -= value;
01362       }
01363       else {
01364         tmp_b_[currentRHS_][row-localOffset_] -= value;
01365         tmp_bc_[row-localOffset_] -= value;
01366       }
01367 
01368       if (debugOutput_) {
01369         fprintf(debugFile_,"BC mod, rhs %d  -= %e\n", row, value);
01370         fprintf(debugFile_,"BC, set A(%d,%d)==%e, to 0.0\n",
01371                 row, bcEqns[idx], offDiagCoefs2[j]);
01372       }
01373 
01374       offDiagCoefs2[j] = 0.0;
01375 
01376     }//for offDiagLength2
01377 
01378   }//for localEnd
01379 
01380   return(0);
01381 }
01382 
01383 //==============================================================================
01384 int Aztec_LinSysCore::enforceBlkEssentialBC(int* blkEqn,
01385                                             int* blkOffset,
01386                                             double* alpha,
01387                                             double* gamma,
01388                                             int len)
01389 {
01390 //
01391 //This function must enforce an essential boundary condition on each local
01392 //equation specified by the pair 'blkEqn' and 'blkOffset'. A point-equation
01393 //resides within a block-equation. The blkOffset gives the distance from the
01394 //beginning of the block-equation to the point-equation.
01395 
01396 //The following modification should be made to A and b, for each incoming
01397 //point-equation:
01398 //
01399 //for(each local equation i){
01400 //   for(each column j in row i) {
01401 //      if (i==j) b[i] = gamma/alpha;
01402 //      else b[j] -= (gamma/alpha)*A[j,i];
01403 //   }
01404 //}
01405 //
01406 //all of row and column 'i' in A should be zeroed,
01407 //except for 1.0 on the diagonal.
01408 //
01409    int val_length = 0;
01410    double* val = NULL;
01411    int colInd_length = 0;
01412    int* blkCols = NULL;
01413    int rowNNZs = 0, numCols = 0, err = 0;
01414 
01415    double* val2 = NULL;
01416    int val2Len = val_length;
01417    int* cols2 = NULL;
01418    int cols2Len = colInd_length;
01419 
01420    int localEnd = localBlkOffset_ + numLocalEqnBlks_ - 1;
01421 
01422    for(int i=0; i<len; i++) {
01423 
01424       //if blkEqn[i] isn't local, then the processor that owns it
01425       //should be running this code too. Otherwise there's trouble...
01426 
01427       if ((localBlkOffset_ > blkEqn[i]) || (blkEqn[i] > localEnd)){
01428          continue;
01429       }
01430 
01431       err = getBlockRow(blkEqn[i], val, val_length, blkCols, colInd_length,
01432                         numCols, rowNNZs);
01433       if (err) {
01434          fei::console_out() << "Aztec_LSC: ERROR in getBlockRow" << FEI_ENDL;
01435          return(-1);
01436       }
01437 
01438       //now let's do the BC modification to this block-row.
01439 
01440       err = blkRowEssBCMod(blkEqn[i], blkOffset[i], val, blkCols, numCols,
01441                              rowNNZs, alpha[i], gamma[i]);
01442 
01443       blkA_ptr_->putBlockRow(blkEqn[i], val, blkCols, numCols);
01444 
01445       //now let's take advantage of the symmetry of element-contributions to
01446       //do the column-wise part of the BC modification. Since the structure of
01447       //the matrix arising from element contributions is symmetric, we know that
01448       //the column indices in row 'blkEqn' correspond to rows which have entries
01449       //in column 'blkEqn'. So we can modify the column in those rows rather
01450       //than searching all rows in the matrix looking for that column.
01451 
01452       //so let's loop over the block-columns and do the column-wise essBC mod
01453       //to those rows.
01454       
01455       for(int j=0; j<numCols; j++) {
01456 
01457          int col_row = blkCols[j];
01458 
01459          //if this column-index doesn't correspond to a local row, skip it.
01460          if ((localOffset_ > col_row) || (col_row > localEnd)) continue;
01461 
01462          //if this column is the diagonal column, skip it.
01463          if (col_row == blkEqn[i]) continue;
01464 
01465          int thisNNZ = 0;
01466          int thisNumBlks = 0;
01467          err = getBlockRow(col_row, val2, val2Len, cols2, cols2Len,
01468                            thisNumBlks, thisNNZ);
01469 
01470          err = blkColEssBCMod(col_row, blkEqn[i], blkOffset[i], val2, cols2,
01471                               thisNumBlks, thisNNZ, alpha[i], gamma[i]);
01472 
01473          blkA_ptr_->putBlockRow(col_row, val2, cols2, thisNumBlks);
01474 
01475       }// end for(j<rowLength) loop
01476    }
01477 
01478    delete [] val;
01479    delete [] blkCols;
01480    delete [] val2;
01481    delete [] cols2;
01482    return(0);
01483 }
01484 
01485 //==============================================================================
01486 int Aztec_LinSysCore::blkRowEssBCMod(int blkEqn, int blkOffset, double* val,
01487                                        int* blkCols, int numCols, int numPtNNZ,
01488                                        double alpha, double gamma)
01489 {
01490 //
01491 //This function performs an essential boundary-condition modification for a
01492 //single block-row.
01493 //
01494    //we need to know which point-row this block-row corresponds to, so
01495    //we can do stuff to the rhs vector.
01496          
01497    int pointRow = blockRowToPointRow(blkEqn);
01498 
01499    int offset = 0;
01500 
01501    for(int j=0; j<numCols; j++) {
01502 
01503       int err, ptRows = 0, ptCols = 0;
01504       err = blkA_ptr_->getBlockSize(blkEqn, blkCols[j], ptRows, ptCols);
01505 
01506       if (err) {
01507          fei::console_out() << "Aztec_LSC::blkRowEssBCMod: error in getBlockSize" << FEI_ENDL;
01508          return(1);
01509       }
01510 
01511       if (blkCols[j] == blkEqn) {
01512          //this is the diagonal block, so we need to diagonalize the
01513          //'blkOffset'th point-row and point-column, leaving a 1.0 on the
01514          //diagonal.
01515          double bc_term = gamma/alpha;
01516          double rhs_term = bc_term;
01517          if (explicitDirichletBCs_) rhs_term = 0.0;
01518 
01519          int thisOffset = offset;
01520 
01521          for(int jj=0; jj<ptCols; jj++) {
01522             if (jj==blkOffset) {
01523                //if this is the point-column of interest, we move the 
01524                //contents over to the rhs in the appropriate way, except
01525                //for the entry on the row-of-interest, which we just
01526                //set to 1.0.
01527 
01528                for(int row=0; row<ptRows; row++) {
01529                   if (row==blkOffset) {
01530                      if (rhsLoaded_) {
01531                         (*b_ptr_)[pointRow+row] = rhs_term;
01532                         (*bc_)[pointRow+row] = bc_term;
01533                      }
01534                      else {
01535                         tmp_b_[currentRHS_][pointRow+row-localOffset_] = rhs_term;
01536                         tmp_bc_[pointRow+row-localOffset_] = bc_term;
01537                      }
01538                      val[thisOffset+row] = 1.0;
01539                   }
01540                   else {
01541                      if (rhsLoaded_) {
01542                         (*b_ptr_)[pointRow+row] -= val[thisOffset+row]
01543                                                 * bc_term;
01544                         (*bc_)[pointRow+row] -= val[thisOffset+row]
01545                                                 * bc_term;
01546                      }
01547                      else {
01548                         tmp_b_[currentRHS_][pointRow+row-localOffset_] -= val[thisOffset+row]*bc_term;
01549                         tmp_bc_[pointRow+row-localOffset_] -= val[thisOffset+row]*bc_term;
01550                      }
01551                      val[thisOffset+row] = 0.0;
01552                   }
01553                }
01554             }
01555             else {
01556                val[thisOffset+blkOffset] = 0.0;
01557             }
01558 
01559             thisOffset += ptRows;
01560          }
01561       }
01562       else {
01563          //this isn't the diagonal block, so we just want to zero the
01564          //whole 'blkOffset'th point-row in this block.
01565          int thisOffset = offset + blkOffset;
01566          for(int ii=0; ii<ptCols; ii++) {
01567             val[thisOffset] = 0.0;
01568             thisOffset += ptRows;
01569          }
01570       }
01571 
01572       offset += ptRows*ptCols;
01573    }
01574 
01575    return(0);
01576 }
01577 
01578 //==============================================================================
01579 int Aztec_LinSysCore::blkColEssBCMod(int blkRow, int blkEqn, int blkOffset,
01580                                      double* val, int* blkCols, int numCols,
01581                                      int numPtNNZ, double alpha, double gamma)
01582 {
01583 //
01584 //This function does the column-wise modification for an Essential BC, to 
01585 //block column 'blkEqn', to the point-equation 'blkOffset' equations into the
01586 //block, for the block-row contained in val,blkCols.
01587 //
01588 //NOTE: blkEqn is a 0-based equation number.
01589 
01590    int thisPtRow = blockRowToPointRow(blkRow);
01591    int thisRowSize = 0, thisColSize = 0;
01592 
01593    //look through the row to find the non-zero blk in position
01594    //blkEqn and make the appropriate modification.
01595    int err, offset = 0;
01596    for(int j=0; j<numCols; j++) {
01597       err = blkA_ptr_->getBlockSize(blkRow, blkCols[j],
01598                                     thisRowSize, thisColSize);
01599 
01600       if (err) {
01601          fei::console_out() << "Aztec_LinSysCore::blkColEssBCMod: ERROR in getBlockSize" << FEI_ENDL;
01602          return(1);
01603       }
01604 
01605       if (blkCols[j] == blkEqn) {
01606          double bc_term = gamma/alpha;
01607 
01608          int thisOffset = offset + blkOffset*thisRowSize;
01609 
01610          for(int row=0; row<thisRowSize; row++) {
01611             if (rhsLoaded_) {
01612                (*b_ptr_)[thisPtRow+row] -= val[thisOffset+row] * bc_term;
01613                (*bc_)[thisPtRow+row] -= val[thisOffset+row] * bc_term;
01614             }
01615             else {
01616                tmp_b_[currentRHS_][thisPtRow+row-localOffset_] -=
01617                                          val[thisOffset+row]*bc_term;
01618                tmp_bc_[thisPtRow+row-localOffset_] -=
01619                                          val[thisOffset+row]*bc_term;
01620             }
01621             val[thisOffset+row] = 0.0;
01622          }
01623 
01624          break;
01625       }
01626 
01627       offset += thisRowSize*thisColSize;
01628    }// end for(j<numCols) loop
01629 
01630    return(0);
01631 }
01632 
01633 //==============================================================================
01634 int Aztec_LinSysCore::blockRowToPointRow(int blkRow) {
01635 //
01636 //This function returns a (global 0-based) point-equation which corresponds to
01637 //the first point-equation in block-row 'blkRow'.
01638 //
01639    int localBlkRow = blkRow - localBlkOffset_;
01640 
01641    if (localBlkRow < 0 || localBlkRow > numLocalEqnBlks_) return(-1);
01642  
01643    int localPtRow = 0;
01644    for(int i=0; i<localBlkRow; i++) {
01645       localPtRow += localBlockSizes_[i];
01646    }
01647 
01648    int pointRow = localPtRow + localOffset_;
01649    return(pointRow);
01650 }
01651 
01652 //==============================================================================
01653 int Aztec_LinSysCore::getBlockRow(int blkRow, double*& val, int& valLen,
01654                                   int*& blkColInds, int& blkColIndLen,
01655                                   int& numNzBlks, int& numNNZ) {
01656 //
01657 //This function gets a block row from the VBR matrix. The val array and the
01658 //blkColInds array are of lengths valLen and blkColIndLen, respectively.
01659 //On output, val and blkColInds are lengthened if necessary, with the new
01660 //lengths updated in valLen and blkColIndLen. The actual number of nonzero
01661 //blocks and nonzero point-entries are returned in numNzBlks and numNNZ.
01662 //
01663    numNNZ = 0;
01664    int err = blkA_ptr_->getNumNonzerosPerRow(blkRow, numNNZ);
01665    if (err) {
01666       fei::console_out() << "Aztec_LSC::getBlockRow: ERROR in getNumNonzerosPerRow" << FEI_ENDL;
01667       return(1);
01668    }
01669 
01670    numNzBlks = 0;
01671    err = blkA_ptr_->getNumBlocksPerRow(blkRow, numNzBlks);
01672    if (err) {
01673       fei::console_out() << "Aztec_LSC::getBlockRow: ERROR in getNumBlocksPerRow" << FEI_ENDL;
01674       return(1);
01675    }
01676 
01677    if (numNNZ > valLen) {
01678       double* newVals = new double[numNNZ];
01679       delete [] val;
01680       val = newVals;
01681       valLen = numNNZ;
01682    }
01683 
01684    if (numNzBlks > blkColIndLen) {
01685       int* newCols = new int[numNzBlks];
01686       delete [] blkColInds;
01687       blkColInds = newCols;
01688       blkColIndLen = numNzBlks;
01689    }
01690 
01691    err = blkA_ptr_->getBlockRow(blkRow, val, blkColInds, numNzBlks);
01692 
01693    return(err);
01694 }
01695 
01696 //==============================================================================
01697 int Aztec_LinSysCore::enforceRemoteEssBCs(int numEqns, int* globalEqns,
01698                                           int** colIndices, int* colIndLen,
01699                                           double** coefs)
01700 {
01701   bcsLoaded_ = true;
01702 
01703   //writeA("preRemBC");
01704 
01705   if (debugOutput_) {
01706     for(int i=0; i<numEqns; ++i) {
01707       fprintf(debugFile_,"remBC row %d, (cols,coefs): ", globalEqns[i]);
01708       for(int j=0; j<colIndLen[i]; ++j) {
01709         fprintf(debugFile_, "(%d,%e) ",colIndices[i][j], coefs[i][j]);
01710       }
01711       fprintf(debugFile_,"\n");
01712     }
01713     fflush(debugFile_);
01714   }
01715 
01716    if (blockMatrix_) {
01717       int* blkEqns = new int[numEqns];
01718       int* blkOffsets = new int[numEqns];
01719 
01720       getBlkEqnsAndOffsets(globalEqns, blkEqns, blkOffsets, numEqns);
01721 
01722       int** blkCols = new int*[numEqns];
01723       int** blkColOffsets = new int*[numEqns];
01724 
01725       for(int i=0; i<numEqns; i++) {
01726          blkCols[i] = new int[colIndLen[i]];
01727          blkColOffsets[i] = new int[colIndLen[i]];
01728 
01729          getBlkEqnsAndOffsets(colIndices[i], blkCols[i], blkColOffsets[i],
01730                               colIndLen[i]);
01731       }
01732 
01733       enforceBlkRemoteEssBCs(numEqns, blkEqns, blkCols, blkColOffsets,
01734                              colIndLen, coefs);
01735 
01736       delete [] blkEqns;
01737       delete [] blkOffsets;
01738 
01739       for(int j=0; j<numEqns; j++) {
01740          delete [] blkCols[j];
01741          delete [] blkColOffsets[j];
01742       }
01743       delete [] blkCols;
01744       delete [] blkColOffsets;
01745 
01746       return(0);
01747    }
01748 
01749    int localEnd = localOffset_ + numLocalEqns_ - 1;
01750 
01751    for(int i=0; i<numEqns; i++) {
01752      int globalEqn_i = globalEqns[i];
01753 
01754      if ((localOffset_ > globalEqn_i) || (globalEqn_i > localEnd)){
01755        continue;
01756      }
01757 
01758      int rowLen = 0;
01759      int* AcolInds = NULL;
01760      double* Acoefs = NULL;
01761 
01762      A_ptr_->getOffDiagRowPointers(globalEqn_i, AcolInds, Acoefs, rowLen);
01763 
01764      for(int j=0; j<colIndLen[i]; j++) {
01765        for(int k=0; k<rowLen; k++) {
01766          if (A_ptr_->getAztec_Map()->getTransformedEqn(AcolInds[k]) == colIndices[i][j]) {
01767            double value = Acoefs[k]*coefs[i][j];
01768 
01769            double old_rhs_val = 0.0;
01770            if (rhsLoaded_) {
01771              old_rhs_val = (*b_ptr_)[globalEqn_i];
01772              (*b_ptr_)[globalEqn_i] -= value;
01773              (*bc_)[globalEqn_i] -= value;
01774            }
01775            else {
01776              old_rhs_val = tmp_b_[currentRHS_][globalEqn_i -localOffset_];
01777              tmp_b_[currentRHS_][globalEqn_i -localOffset_] -= value;
01778              tmp_bc_[globalEqn_i -localOffset_] -= value;
01779            }
01780 
01781            if (debugOutput_) {
01782              fprintf(debugFile_,"remBC mod, rhs %d (%e) -= %e\n",
01783                      globalEqn_i, old_rhs_val, value);
01784              fprintf(debugFile_,"remBC, set A(%d,%d)==%e, to 0.0\n",
01785                      globalEqn_i, AcolInds[k], Acoefs[k]);
01786            }
01787 
01788            Acoefs[k] = 0.0;
01789          }
01790        }
01791      }
01792 
01793    }
01794 
01795    return(0);
01796 }
01797 
01798 //==============================================================================
01799 int Aztec_LinSysCore::enforceBlkRemoteEssBCs(int numEqns, int* blkEqns,
01800                                int** blkColInds, int** blkColOffsets,
01801                                int* blkColLens, double** remEssBCCoefs) {
01802    int val_length = 0;
01803    double* val = NULL;
01804    int colInd_length = 0;
01805    int* blkCols = NULL;
01806    int rowNNZs = 0, numCols = 0, err = 0;
01807 
01808    int localEnd = localBlkOffset_ + numLocalEqnBlks_ - 1;
01809 
01810    for(int i=0; i<numEqns; i++) {
01811       if ((localBlkOffset_ > blkEqns[i]) || (blkEqns[i] > localEnd)){
01812          continue;
01813       }
01814 
01815       err = getBlockRow(blkEqns[i], val, val_length, blkCols, colInd_length,
01816                         numCols, rowNNZs);
01817       if (err) {
01818          fei::console_out() << "Aztec_LSC:enforceBlkRemoteEssBC ERROR in getBlockRow" << FEI_ENDL;
01819          return(-1);
01820       }
01821 
01822       //we need to know which point-row this block-row corresponds to, so
01823       //we can do stuff to the rhs vector.
01824 
01825       int pointRow = blockRowToPointRow(blkEqns[i]);
01826 
01827       int offset = 0;
01828 
01829       for(int j=0; j<numCols; j++) {
01830          int ptRows = 0, ptCols = 0;
01831          err = blkA_ptr_->getBlockSize(blkEqns[i], blkCols[j],
01832                                        ptRows, ptCols);
01833          if (err) {
01834             fei::console_out() << "Aztec_LSC::enforceBlkRemoteEssBCs: error in getBlockSize"
01835                 << FEI_ENDL;
01836             return(-1);
01837          }
01838 
01839          for(int k=0; k<blkColLens[i]; k++) {
01840             if (blkColInds[i][k] == blkCols[j]) {
01841                int thisOffset = offset + blkColOffsets[i][k] * ptRows;
01842                double rhsTerm = remEssBCCoefs[i][k];
01843 
01844                double* bvec = &(tmp_b_[currentRHS_][pointRow-localOffset_]);
01845                double* bcvec = &(tmp_bc_[pointRow-localOffset_]);
01846                for(int row=0; row<ptRows; row++) {
01847                   double& coef = val[thisOffset+row];
01848                   bvec[row] -= coef*rhsTerm;
01849                   bcvec[row] -= coef*rhsTerm;
01850                   coef = 0.0;
01851                }
01852 
01853                blkA_ptr_->putBlockRow(blkEqns[i], &val[offset],
01854                                       &blkCols[j], 1);
01855             }
01856          }
01857 
01858          offset += ptRows*ptCols;
01859       }
01860    }
01861 
01862    delete [] val;
01863    delete [] blkCols;
01864    return(0);
01865 }
01866 
01867 //==============================================================================
01868 int Aztec_LinSysCore::getMatrixPtr(Data& data)
01869 {
01870    if (!matrixLoaded_) matrixLoadComplete();
01871 
01872    if (blockMatrix_) {
01873       data.setTypeName("AztecDVBR_Matrix");
01874       data.setDataPtr((void*)blkA_ptr_);
01875    }
01876    else {
01877       data.setTypeName("AztecDMSR_Matrix");
01878       data.setDataPtr((void*)A_ptr_);
01879    }
01880    return(0);
01881 }
01882 
01883 //==============================================================================
01884 int Aztec_LinSysCore::copyInMatrix(double scalar, const Data& data) {
01885 //
01886 //Overwrites the current internal matrix with a scaled copy of the
01887 //input argument.
01888 //
01889    if (blockMatrix_) {
01890       if (strcmp("AztecDVBR_Matrix", data.getTypeName()))
01891          messageAbort("copyInMatrix: data type string not 'AztecDVBR_Matrix'.");
01892 
01893       AztecDVBR_Matrix* source = (AztecDVBR_Matrix*)data.getDataPtr();
01894 
01895       if (blkA_ != NULL) delete blkA_;
01896       blkA_ = new AztecDVBR_Matrix(*source);
01897       blkA_ptr_ = blkA_;
01898 
01899       VBRmatPlusScaledMat(blkA_ptr_, scalar, source);
01900 
01901       azA_ = blkA_ptr_->getAZ_MATRIX_Ptr();
01902    }
01903    else {
01904       if (strcmp("AztecDMSR_Matrix", data.getTypeName()))
01905          messageAbort("copyInMatrix: data type string not 'AztecDMSR_Matrix'.");
01906 
01907       AztecDMSR_Matrix* source = (AztecDMSR_Matrix*)data.getDataPtr();
01908       A_ptr_->copyStructure(*source);
01909 
01910       MSRmatPlusScaledMat(A_ptr_, scalar, source);
01911 
01912       azA_ = A_ptr_->getAZ_MATRIX_PTR();
01913    }
01914 
01915    needNewPreconditioner_ = true;
01916    return(0);
01917 }
01918 
01919 //==============================================================================
01920 int Aztec_LinSysCore::copyOutMatrix(double scalar, Data& data) {
01921 //
01922 //Passes out a scaled copy of the current internal matrix.
01923 //
01924 
01925    if (!matrixLoaded_) matrixLoadComplete();
01926 
01927    if (blockMatrix_) {
01928       AztecDVBR_Matrix* outmat = new AztecDVBR_Matrix(*blkA_ptr_);
01929 
01930       //the copy-constructor above took all structural info from blkA_ptr_,
01931       //but not the data coefficients.
01932 
01933       VBRmatPlusScaledMat(outmat, scalar, blkA_ptr_);
01934 
01935       data.setTypeName("AztecDVBR_Matrix");
01936       data.setDataPtr((void*)outmat);
01937    }
01938    else {
01939       AztecDMSR_Matrix* outmat = new AztecDMSR_Matrix(*A_ptr_);
01940 
01941       outmat->scale(scalar);
01942 
01943       data.setTypeName("AztecDMSR_Matrix");
01944       data.setDataPtr((void*)outmat);
01945    }
01946    return(0);
01947 }
01948 
01949 //==============================================================================
01950 int Aztec_LinSysCore::sumInMatrix(double scalar, const Data& data) {
01951 
01952    if (!matrixLoaded_) matrixLoadComplete();
01953 
01954    if (blockMatrix_) {
01955      if (strcmp("AztecDVBR_Matrix", data.getTypeName())) {
01956        fei::console_out() << "Aztec_LinSysCore::sumInMatrix ERROR, incoming type-string: "
01957             << data.getTypeName() << ", expected AztecDVBR_Matrix." << FEI_ENDL;
01958        messageAbort("Aborting.");
01959      }
01960      AztecDVBR_Matrix* source = (AztecDVBR_Matrix*)data.getDataPtr();
01961 
01962      VBRmatPlusScaledMat(blkA_ptr_, scalar, source);
01963    }
01964    else {
01965       if (strcmp("AztecDMSR_Matrix", data.getTypeName()))
01966          messageAbort("sumInMatrix: data type string not 'AztecDMSR_Matrix'.");
01967 
01968       AztecDMSR_Matrix* source = (AztecDMSR_Matrix*)data.getDataPtr();
01969 
01970       MSRmatPlusScaledMat(A_ptr_, scalar, source);
01971    }
01972 
01973    needNewPreconditioner_ = true;
01974    return(0);
01975 }
01976 
01977 //==============================================================================
01978 int Aztec_LinSysCore::getRHSVectorPtr(Data& data) {
01979 
01980    if (!matrixLoaded_) matrixLoadComplete();
01981 
01982    data.setTypeName("Aztec_LSVector");
01983    data.setDataPtr((void*)b_ptr_);
01984    return(0);
01985 }
01986 
01987 //==============================================================================
01988 int Aztec_LinSysCore::copyInRHSVector(double scalar, const Data& data) {
01989 
01990    if (!rhsLoaded_) matrixLoadComplete();
01991 
01992    if (strcmp("Aztec_LSVector", data.getTypeName()))
01993       messageAbort("copyInRHSVector: data's type string not 'Aztec_LSVector'.");
01994 
01995    Aztec_LSVector* sourcevec = (Aztec_LSVector*)data.getDataPtr();
01996 
01997    *b_ptr_ = *sourcevec;
01998 
01999    if (scalar != 1.0) b_ptr_->scale(scalar);
02000    return(0);
02001 }
02002 
02003 //==============================================================================
02004 int Aztec_LinSysCore::copyOutRHSVector(double scalar, Data& data) {
02005 
02006    if (!rhsLoaded_) matrixLoadComplete();
02007 
02008    Aztec_LSVector* outvec = new Aztec_LSVector(*b_ptr_);
02009 
02010    outvec->put(0.0);
02011 
02012    outvec->addVec(scalar, *b_ptr_);
02013 
02014    data.setTypeName("Aztec_LSVector");
02015    data.setDataPtr((void*)outvec);
02016    return(0);
02017 }
02018 
02019 //==============================================================================
02020 int Aztec_LinSysCore::sumInRHSVector(double scalar, const Data& data) {
02021 
02022    if (!rhsLoaded_) matrixLoadComplete();
02023 
02024    if (strcmp("Aztec_LSVector", data.getTypeName()))
02025       messageAbort("sumInRHSVector: data's type string not 'Aztec_LSVector'.");
02026 
02027    Aztec_LSVector* source = (Aztec_LSVector*)data.getDataPtr();
02028 
02029    b_ptr_->addVec(scalar, *source);
02030    return(0);
02031 }
02032 
02033 //==============================================================================
02034 int Aztec_LinSysCore::destroyMatrixData(Data& data) {
02035 
02036    if (blockMatrix_) {
02037       if (strcmp("AztecDVBR_Matrix", data.getTypeName()))
02038          messageAbort("destroyMatrixData: data doesn't contain a AztecDVBR_Matrix.");
02039 
02040       AztecDVBR_Matrix* mat = (AztecDVBR_Matrix*)data.getDataPtr();
02041       delete mat;
02042    }
02043    else {
02044       if (strcmp("AztecDMSR_Matrix", data.getTypeName()))
02045          messageAbort("destroyMatrixData: data doesn't contain a AztecDMSR_Matrix.");
02046 
02047       AztecDMSR_Matrix* mat = (AztecDMSR_Matrix*)data.getDataPtr();
02048       delete mat;
02049    }
02050    return(0);
02051 }
02052 
02053 //==============================================================================
02054 int Aztec_LinSysCore::destroyVectorData(Data& data) {
02055 
02056    if (strcmp("Aztec_LSVector", data.getTypeName()))
02057       messageAbort("destroyVectorData: data doesn't contain a Aztec_LSVector.");
02058 
02059    Aztec_LSVector* vec = (Aztec_LSVector*)data.getDataPtr();
02060    delete vec;
02061    return(0);
02062 }
02063 
02064 //==============================================================================
02065 int Aztec_LinSysCore::setNumRHSVectors(int numRHSs, const int* rhsIDs) {
02066 
02067    if (numRHSs < 0)
02068       messageAbort("setNumRHSVectors: numRHSs < 0.");
02069 
02070    if (numRHSs == 0) return(0);
02071 
02072    delete [] rhsIDs_;
02073    numRHSs_ = numRHSs;
02074    rhsIDs_ = new int[numRHSs_];
02075 
02076    for(int i=0; i<numRHSs_; i++) rhsIDs_[i] = rhsIDs[i];
02077    return(0);
02078 }
02079 
02080 //==============================================================================
02081 int Aztec_LinSysCore::setRHSID(int rhsID) {
02082 
02083    for(int i=0; i<numRHSs_; i++){
02084       if (rhsIDs_[i] == rhsID){
02085          currentRHS_ = i;
02086          if (rhsLoaded_) b_ptr_ = b_[currentRHS_];
02087          return(0);
02088       }
02089    }
02090 
02091    messageAbort("setRHSID: rhsID not found.");
02092    return(-1);
02093 }
02094 
02095 //==============================================================================
02096 int Aztec_LinSysCore::putInitialGuess(const int* eqnNumbers,
02097                                       const double* values,
02098                                       int len) {
02099 //
02100 //This function scatters (puts) values into the linear-system's soln vector.
02101 //
02102 // num is how many values are being passed,
02103 // indices holds the global 'row-numbers' into which the values go,
02104 // and values holds the actual coefficients to be scattered.
02105 //
02106 
02107    int localEnd = localOffset_ + numLocalEqns_ -1;
02108    if (matrixLoaded_) {
02109       for(int i=0; i<len; i++){
02110          if ((localOffset_ > eqnNumbers[i]) || (eqnNumbers[i] > localEnd))
02111             continue;
02112 
02113          (*x_)[eqnNumbers[i]] = values[i];
02114       }
02115    }
02116    else {
02117       tmp_x_touched_ = true;
02118       for(int i=0; i<len; i++){
02119          if ((localOffset_ > eqnNumbers[i]) || (eqnNumbers[i] > localEnd))
02120             continue;
02121          tmp_x_[eqnNumbers[i]-localOffset_] = values[i];
02122       }
02123    }
02124    return(0);
02125 }
02126 
02127 //==============================================================================
02128 int Aztec_LinSysCore::getSolution(double* answers, int len) {
02129 //
02130 //The caller must allocate the memory for 'answers',
02131 //and len must be set to the right value -- i.e., len should equal
02132 //numLocalEqns_.
02133 //
02134    if (len != numLocalEqns_)
02135       messageAbort("getSolution: len != numLocalEqns_.");
02136 
02137    for(int i=0; i<numLocalEqns_; i++) {
02138       answers[i] = (*x_)[localOffset_ + i];
02139    }
02140    return(0);
02141 }
02142 
02143 //==============================================================================
02144 int Aztec_LinSysCore::getSolnEntry(int eqnNumber, double& answer) {
02145 //
02146 //This function returns a single solution entry, the coefficient for
02147 //equation number eqnNumber.
02148 //
02149    int localEnd = localOffset_ + numLocalEqns_ -1;
02150    if ((localOffset_ > eqnNumber) || (eqnNumber > localEnd))
02151       messageAbort("getSolnEntry: eqnNumber out of range.");
02152 
02153    answer = (*x_)[eqnNumber];
02154    return(0);
02155 }
02156 
02157 //==============================================================================
02158 int Aztec_LinSysCore::formResidual(double* values, int len)
02159 {
02160    if (len != numLocalEqns_) {
02161       messageAbort("formResidual: len != numLocalEqns_.");
02162    }
02163 
02164    if (!matrixLoaded_ || !rhsLoaded_) {
02165       matrixLoadComplete();
02166    }
02167 
02168    //after the solve x_ and b_ptr_ were 'inv_order'd back into user-ordering.
02169    //Before we can do this residual calculation, we have to put them into Aztec
02170    //ordering again.
02171 
02172    int* update_index = NULL;
02173    if (blockMatrix_) {
02174       update_index = blkA_ptr_->getUpdate_index();
02175    }
02176    else {
02177       update_index = A_ptr_->getAztec_Map()->update_index;
02178    }
02179 
02180    AZ_reorder_vec((double*)(x_->startPointer()), azA_->data_org, update_index,
02181                   azA_->rpntr);
02182 
02183    AZ_reorder_vec((double*)(b_ptr_->startPointer()), azA_->data_org,
02184                   update_index, azA_->rpntr);
02185 
02186    Aztec_LSVector* r = new Aztec_LSVector(*x_);
02187 
02188    if (blockMatrix_) blkA_ptr_->matvec(*x_, *r); //form r = A*x
02189    else A_ptr_->matvec(*x_, *r);
02190 
02191    r->addVec(-1.0, *b_ptr_); //form r = A*x - b
02192 
02193    r->scale(-1.0); //form r = b - A*x (not really necessary, but...)
02194 
02195    //now let's get the residual r into user ordering...
02196 
02197    Aztec_LSVector* rtmp = new Aztec_LSVector(*x_);
02198 
02199    AZ_invorder_vec((double*)(r->startPointer()), azA_->data_org, update_index,
02200                    azA_->rpntr, (double*)rtmp->startPointer());
02201 
02202    *r = *rtmp;
02203 
02204    const double* rptr = r->startPointer();
02205 
02206    for(int i=0; i<numLocalEqns_; i++) {
02207       values[i] = rptr[i];
02208    }
02209 
02210    //finally, let's put x_ and b_ptr_ back into user ordering before we exit...
02211    AZ_invorder_vec((double*)(x_->startPointer()), azA_->data_org, update_index,
02212                    azA_->rpntr, (double*)rtmp->startPointer());
02213 
02214    *x_ = *rtmp;
02215 
02216    AZ_invorder_vec((double*)(b_ptr_->startPointer()), azA_->data_org,
02217                    update_index, azA_->rpntr, (double*)rtmp->startPointer());
02218 
02219    *b_ptr_ = *rtmp;
02220 
02221    delete rtmp;
02222    delete r;
02223    return(0);
02224 }
02225 
02226 //==============================================================================
02227 int Aztec_LinSysCore::selectSolver(const char* name) {
02228 //
02229 // This function takes a string naming the solver and sets the solver choice
02230 // in the options array accordingly.
02231 //
02232    char msg[64];
02233 
02234    if (!strcmp(name, "AZ_gmres")) {
02235       aztec_options_[AZ_solver] = AZ_gmres;
02236       sprintf(msg, "AZ_gmres solver.");
02237    }
02238    else if (!strcmp(name, "AZ_cg")) {
02239       aztec_options_[AZ_solver] = AZ_cg;
02240       sprintf(msg, "AZ_cg solver.");
02241    }
02242    else if (!strcmp(name, "AZ_bicgstab")) {
02243       aztec_options_[AZ_solver] = AZ_bicgstab;
02244       sprintf(msg, "AZ_bicgstab solver.");
02245    }
02246    else if (!strcmp(name, "AZ_cgs")) {
02247       aztec_options_[AZ_solver] = AZ_cgs;
02248       sprintf(msg, "AZ_cgs solver.");
02249    }
02250    else if (!strcmp(name, "AZ_tfqmr")) {
02251       aztec_options_[AZ_solver] = AZ_tfqmr;
02252       sprintf(msg, "AZ_tfqmr solver.");
02253    }
02254    else if (!strcmp(name, "AZ_lu")) {
02255       aztec_options_[AZ_solver] = AZ_lu;
02256       sprintf(msg, "AZ_lu solver.");
02257    }
02258    else {
02259       aztec_options_[AZ_solver] = AZ_gmres;
02260       sprintf(msg, "AZ_gmres default solver.");
02261       if (thisProc_ == 0) {
02262          FEI_COUT << "Aztec_LinSysCore: Warning: requested solver <" << name << "> not recognized, defaulting to AZ_gmres." << FEI_ENDL;
02263       }
02264    }
02265 
02266    debugOutput(msg);
02267 
02268    return(0);
02269 }
02270 
02271 //==============================================================================
02272 int Aztec_LinSysCore::selectPreconditioner(const char* name) {
02273 
02274    char msg[128];
02275    sprintf(msg, "selectPreconditioner(%s)", name);
02276    debugOutput(msg);
02277 
02278    if (!strcmp(name, "AZ_none")) {
02279       aztec_options_[AZ_precond] = AZ_none;
02280       sprintf(msg, " -- selected: AZ_none.");
02281    }
02282    else if (!strcmp(name, "AZ_Jacobi")) {
02283       aztec_options_[AZ_precond] = AZ_Jacobi;
02284       sprintf(msg, " -- selected: AZ_Jacobi.");
02285    }
02286    else if (!strcmp(name, "AZ_Neumann")) {
02287       aztec_options_[AZ_precond] = AZ_Neumann;
02288       sprintf(msg, " -- selected: AZ_Neumann.");
02289    }
02290    else if (!strcmp(name, "AZ_ls")) {
02291       aztec_options_[AZ_precond] = AZ_ls;
02292       sprintf(msg, " -- selected: AZ_ls.");
02293    }
02294    else if (!strcmp(name, "AZ_sym_GS")) {
02295       aztec_options_[AZ_precond] = AZ_sym_GS;
02296       sprintf(msg, " -- selected: AZ_sym_GS.");
02297    }
02298    else if (!strcmp(name, "AZ_dom_decomp")) {
02299       aztec_options_[AZ_precond] = AZ_dom_decomp;
02300       sprintf(msg, " -- selected: AZ_dom_decomp.");
02301    }
02302 //#ifndef NOT_USING_ML
02303 #ifdef USING_ML
02304    else if (!strcmp(name, "ML_Vanek")) {
02305       ML_Vanek_ = true;
02306       aztec_options_[AZ_precond] = AZ_user_precond;
02307       sprintf(msg, " -- selected: AZ_user_precond.");
02308    }
02309 #endif
02310    else {
02311       aztec_options_[AZ_precond] = AZ_none;
02312       sprintf(msg," -- selected: Default, AZ_none.\n");
02313       if (thisProc_ == 0) {
02314          FEI_COUT << "Aztec_LinSysCore: Warning: requested preconditioner <" << name << "> not recognized, defaulting to AZ_none." << FEI_ENDL;
02315       }
02316    }
02317 
02318    debugOutput(msg);
02319 
02320    return(0);
02321 }
02322 
02323 //==============================================================================
02324 void Aztec_LinSysCore::setSubdomainSolve(const char* name) {
02325 
02326    char msg[128];
02327    sprintf(msg, "setSubdomainSolve(%s)", name);
02328    debugOutput(msg);
02329 
02330    if (!strcmp(name, "AZ_lu")) {
02331       aztec_options_[AZ_subdomain_solve] = AZ_lu;
02332       sprintf(msg, " -- selected AZ_lu");
02333    }
02334    else if (!strcmp(name, "AZ_ilu")) {
02335       aztec_options_[AZ_subdomain_solve] = AZ_ilu;
02336       sprintf(msg, " -- selected AZ_ilu");
02337    }
02338    else if (!strcmp(name, "AZ_ilut")) {
02339       aztec_options_[AZ_subdomain_solve] = AZ_ilut;
02340       sprintf(msg, " -- selected AZ_ilut");
02341    }
02342    else if (!strcmp(name, "AZ_rilu")) {
02343       aztec_options_[AZ_subdomain_solve] = AZ_rilu;
02344       sprintf(msg, " -- selected AZ_rilu");
02345    }
02346    else if (!strcmp(name, "AZ_bilu")) {
02347       aztec_options_[AZ_subdomain_solve] = AZ_bilu;
02348       sprintf(msg, " -- selected AZ_bilu");
02349    }
02350    else if (!strcmp(name, "AZ_icc")) {
02351       aztec_options_[AZ_subdomain_solve] = AZ_icc;
02352       sprintf(msg, " -- selected AZ_icc");
02353    }
02354    else {
02355       if (thisProc_ == 0) {
02356          FEI_COUT << "Aztec_LinSysCore: Warning: requested subdomain-solve <" << name << "> not recognized." << FEI_ENDL;
02357       }
02358    }
02359 
02360    debugOutput(msg);
02361 }
02362 
02363 //==============================================================================
02364 void Aztec_LinSysCore::setTypeOverlap(const char* name) {
02365 
02366    char msg[128];
02367    sprintf(msg, "setTypeOverlap(%s)", name);
02368    debugOutput(msg);
02369 
02370    if (!strcmp(name, "AZ_standard")) {
02371       aztec_options_[AZ_type_overlap] = AZ_standard;
02372       sprintf(msg, " -- selected AZ_standard");
02373    }
02374    else if (!strcmp(name, "AZ_ilu")) {
02375       aztec_options_[AZ_type_overlap] = AZ_symmetric;
02376       sprintf(msg, " -- selected AZ_symmetric");
02377    }
02378    else {
02379       if (thisProc_ == 0) {
02380          FEI_COUT << "Aztec_LinSysCore: Warning: requested type-overlap <" << name << "> not recognized." << FEI_ENDL;
02381       }
02382    }
02383 
02384    debugOutput(msg);
02385 }
02386 
02387 //==============================================================================
02388 int Aztec_LinSysCore::writeSystem(const char* name)
02389 {
02390    writeA(name);
02391    return(0);
02392 }
02393 
02394 //==============================================================================
02395 void Aztec_LinSysCore::recordUserParams()
02396 {
02397    checkForParam("AZ_tol", numParams_, paramStrings_,
02398                   aztec_params_[AZ_tol]);
02399 
02400    checkForParam("AZ_drop", numParams_, paramStrings_,
02401                   aztec_params_[AZ_drop]);
02402 
02403    checkForParam("AZ_ilut_fill", numParams_, paramStrings_,
02404                  aztec_params_[AZ_ilut_fill]);
02405 
02406    checkForParam("AZ_omega", numParams_, paramStrings_,
02407                  aztec_params_[AZ_omega]);
02408 
02409    checkForParam("AZ_weights", numParams_, paramStrings_,
02410                  aztec_params_[AZ_weights]);
02411 }
02412 
02413 //==============================================================================
02414 void Aztec_LinSysCore::recordUserOptions()
02415 {
02416 //
02417 //Private function, to be called from launchSolver.
02418 //
02419    const char* param = NULL;
02420 
02421    param = snl_fei::getParamValue("AZ_solver", numParams_, paramStrings_);
02422    if (param != NULL){
02423       selectSolver(param);
02424    }
02425  
02426    param = snl_fei::getParamValue("AZ_precond", numParams_, paramStrings_);
02427    if (param != NULL){
02428       selectPreconditioner(param);
02429    }
02430 
02431    param = snl_fei::getParamValue("AZ_subdomain_solve",
02432                                          numParams_, paramStrings_);
02433    if (param != NULL){
02434       setSubdomainSolve(param);
02435    }
02436 
02437    param = snl_fei::getParamValue("AZ_scaling", numParams_, paramStrings_);
02438    if (param != NULL){
02439       setScalingOption(param);
02440    }
02441 
02442    param = snl_fei::getParamValue("AZ_conv", numParams_, paramStrings_);
02443    if (param != NULL){
02444       setConvTest(param);
02445    }
02446 
02447    param = snl_fei::getParamValue("AZ_pre_calc",
02448                                          numParams_, paramStrings_);
02449    if (param != NULL){
02450       setPreCalc(param);
02451    }
02452 
02453    param = snl_fei::getParamValue("AZ_overlap", numParams_, paramStrings_);
02454    if (param != NULL){
02455       setOverlap(param);
02456    }
02457 
02458    param = snl_fei::getParamValue("AZ_type_overlap",
02459                                          numParams_, paramStrings_);
02460    if (param != NULL){
02461       setTypeOverlap(param);
02462    }
02463 
02464    param = snl_fei::getParamValue("AZ_orthog", numParams_, paramStrings_);
02465    if (param != NULL){
02466       setOrthog(param);
02467    }
02468 
02469    param = snl_fei::getParamValue("AZ_aux_vec", numParams_, paramStrings_);
02470    if (param != NULL){
02471       setAuxVec(param);
02472    }
02473 
02474    param = snl_fei::getParamValue("AZ_output", numParams_, paramStrings_);
02475    if (param != NULL){
02476       setAZ_output(param);
02477    }
02478    else aztec_options_[AZ_output] = outputLevel_;
02479 
02480    checkForOption("AZ_poly_ord", numParams_, paramStrings_,
02481                   aztec_options_[AZ_poly_ord]);
02482 
02483    checkForOption("AZ_kspace", numParams_, paramStrings_,
02484                   aztec_options_[AZ_kspace]);
02485 
02486    checkForOption("AZ_max_iter", numParams_, paramStrings_,
02487                   aztec_options_[AZ_max_iter]);
02488 
02489    checkForOption("AZ_reorder", numParams_, paramStrings_,
02490                   aztec_options_[AZ_reorder]);
02491 
02492    checkForOption("AZ_graph_fill", numParams_, paramStrings_,
02493                   aztec_options_[AZ_graph_fill]);
02494 
02495    checkForOption("AZ_keep_info", numParams_, paramStrings_,
02496                   aztec_options_[AZ_keep_info]);
02497 }
02498 
02499 //==============================================================================
02500 int Aztec_LinSysCore::writeA(const char* name)
02501 {
02502   if (name == NULL) {
02503     return(-1);
02504   }
02505 
02506   FEI_OSTRINGSTREAM osstr;
02507 
02508   if (debugPath_ != NULL) {
02509     osstr << debugPath_;
02510   }
02511   else osstr << ".";
02512 
02513   osstr << "/A_" << name << ".mtx";
02514 
02515   if (blockMatrix_) osstr << ".vbr";
02516 
02517   std::string str = osstr.str();
02518   const char* matname = str.c_str();
02519 
02520   if (blockMatrix_) {
02521     blkA_ptr_->writeToFile(matname);
02522   }
02523   else {
02524     A_ptr_->writeToFile(matname);
02525   }
02526 
02527   return(0);
02528 }
02529 
02530 //==============================================================================
02531 int Aztec_LinSysCore::writeVec(Aztec_LSVector* v, const char* name)
02532 {
02533   if (name == NULL || v == NULL) {
02534     return(-1);
02535   }
02536 
02537   FEI_OSTRINGSTREAM osstr;
02538 
02539   if (debugPath_ != NULL) {
02540     osstr << debugPath_;
02541   }
02542   else osstr << ".";
02543 
02544   osstr << "/" << name << ".vec";
02545 
02546   std::string str = osstr.str();
02547 
02548   const char* vecname = str.c_str();
02549 
02550   v->writeToFile(vecname);
02551 
02552   return(0);
02553 }
02554 
02555 //==============================================================================
02556 int Aztec_LinSysCore::modifyRHSforBCs()
02557 {
02558   for(int i=0; i<numLocalEqns_; i++) {
02559     (*b_ptr_)[i+localOffset_] += tmp_bc_[i];
02560   }
02561 
02562   if (explicitDirichletBCs_) {
02563     for(int j=0; j<numEssBCs_; j++) {
02564       int index = essBCindices_[j];
02565       (*b_ptr_)[index] = 0.0;
02566     }
02567   }
02568   else {
02569     for(int j=0; j<numEssBCs_; j++) {
02570       int index = essBCindices_[j];
02571       (*b_ptr_)[index] = tmp_bc_[index-localOffset_];
02572     }
02573   }
02574 
02575   return(0);
02576 }
02577 
02578 //==============================================================================
02579 int Aztec_LinSysCore::explicitlySetDirichletBCs()
02580 {
02581   for(int j=0; j<numEssBCs_; j++) {
02582     int index = essBCindices_[j];
02583     if (rhsLoaded_) {
02584       (*b_ptr_)[index] = (*bc_)[index];
02585       (*x_)[index] = (*bc_)[index];
02586     }
02587     else {
02588       (*b_ptr_)[index] = tmp_bc_[index-localOffset_];
02589       (*x_)[index] = tmp_bc_[index-localOffset_];
02590     }
02591   }
02592   return(0);
02593 }
02594 
02595 //==============================================================================
02596 int Aztec_LinSysCore::launchSolver(int& solveStatus, int& iterations) {
02597 //
02598 //This function does any last-second setup required for the
02599 //linear solver, then goes ahead and launches the solver to get
02600 //the solution vector.
02601 //Also, if possible, the number of iterations that were performed
02602 //is stored in the iterations_ variable.
02603 //
02604 
02605    unsigned counter = 0;
02606    std::map<std::string,unsigned>::iterator 
02607      iter = named_solve_counter_.find(name_);
02608    if (iter == named_solve_counter_.end()) {
02609      fei::console_out() << "fei: Aztec_LinSysCore::LaunchSolver: internal error."
02610       << FEI_ENDL;
02611    }
02612    else {
02613      counter = iter->second++;
02614    }
02615 
02616    if (debugOutput_ && outputLevel_ > 1) {
02617       FEI_OSTRINGSTREAM osstr;
02618       osstr << name_ << "_Aztec.np"<<numProcs_<<".slv"<<counter;
02619       std::string str = osstr.str();
02620 
02621       writeA(str.c_str());
02622 
02623       FEI_OSTRINGSTREAM x0_osstr;
02624       x0_osstr << "x0_" << str;
02625       std::string x0_str = x0_osstr.str();
02626 
02627       writeVec(x_, x0_str.c_str());
02628 
02629       FEI_OSTRINGSTREAM b_osstr;
02630       b_osstr << "b_" << str;
02631       std::string b_str = b_osstr.str();
02632 
02633       writeVec(b_ptr_, b_str.c_str());
02634    }
02635 
02636    if (needNewPreconditioner_) {
02637      if (precondCreated_) {
02638        AZ_precond_destroy(&azP_);
02639      }
02640 
02641 //#ifndef NOT_USING_ML
02642 #ifdef USING_ML
02643      ML* ml = NULL;
02644 
02645      if (ML_Vanek_) {
02646        int numFineSweeps = 2;
02647        int numCoarseSweeps = 2;
02648        double omega = 0.67;
02649 
02650        initialize_ML(azA_, &azP_, numLevels_,
02651                      numFineSweeps, numCoarseSweeps, omega,
02652                      map_->getProcConfig(), &ml);
02653      }
02654      else {
02655        //set the preconditioning matrix Pmat to point to azA_ (Amat).
02656        azA_->data_org[AZ_name] = 0;
02657        azP_ = AZ_precond_create(azA_, AZ_precondition, NULL);
02658      }
02659 #else
02660      azA_->data_org[AZ_name] = 0;
02661      azP_ = AZ_precond_create(azA_, AZ_precondition, NULL);
02662 #endif
02663      needNewPreconditioner_ = false;
02664      aztec_options_[AZ_pre_calc] = AZ_calc;
02665      aztec_options_[AZ_conv] = AZ_rhs;
02666      aztec_options_[AZ_orthog] = AZ_modified;
02667 
02668      recordUserParams();
02669      recordUserOptions();
02670 
02671      aztec_options_[AZ_keep_info] = 1;
02672    }
02673    else {
02674      aztec_options_[AZ_pre_calc] = AZ_reuse;
02675    }
02676 
02677    precondCreated_ = true;
02678 
02679    int* proc_config = NULL;
02680    int* update_index = NULL;
02681    if (blockMatrix_) {
02682       proc_config = blkMap_->getProcConfig();
02683       update_index = blkA_ptr_->getUpdate_index();
02684    }
02685    else {
02686       proc_config = map_->getProcConfig();
02687       update_index = A_ptr_->getAztec_Map()->update_index;
02688    }
02689 
02690    AZ_reorder_vec((double*)(x_->startPointer()), azA_->data_org, update_index,
02691                   azA_->rpntr);
02692 
02693    AZ_reorder_vec((double*)(b_ptr_->startPointer()), azA_->data_org,
02694                   update_index, azA_->rpntr);
02695 
02696    AZ_iterate((double*)(x_->startPointer()),
02697               (double*)(b_ptr_->startPointer()),
02698               aztec_options_, aztec_params_, aztec_status_,
02699               proc_config, azA_, azP_, azS_);
02700 
02701    iterations = (int)aztec_status_[AZ_its];
02702 
02703    solveStatus = (int)aztec_status_[AZ_why];
02704 
02705    azlsc_solveCounter_++;
02706 
02707    Aztec_LSVector* xtmp = new Aztec_LSVector(*x_);
02708 
02709    //now we need to put x_ back into user-ordering for when we're asked to
02710    //hand out solution entries.
02711 
02712    AZ_invorder_vec((double*)(x_->startPointer()), azA_->data_org, update_index,
02713                    azA_->rpntr, (double*)xtmp->startPointer());
02714 
02715    *x_ = *xtmp;
02716 
02717    //let's put b back into user-ordering too...
02718    AZ_invorder_vec((double*)(b_ptr_->startPointer()), azA_->data_org,
02719                    update_index, azA_->rpntr, (double*)xtmp->startPointer());
02720 
02721    *b_ptr_ = *xtmp;
02722 
02723    delete xtmp;
02724 
02725    if (explicitDirichletBCs_) explicitlySetDirichletBCs();
02726 
02727    if (debugOutput_) {
02728       FEI_OSTRINGSTREAM osstr;
02729       osstr << name_ << "_Aztec.slv"<<counter;
02730       std::string str = osstr.str();
02731 
02732       FEI_OSTRINGSTREAM x_osstr;
02733       x_osstr << "x_" << str;
02734       std::string x_str = x_osstr.str();
02735 
02736       writeVec(x_, x_str.c_str());
02737    }
02738    return(0);
02739 }
02740 
02741 //==============================================================================
02742 void Aztec_LinSysCore::setScalingOption(const char* param) {
02743 
02744    char* msg = new char[128];
02745    for(int i=0; i<128; i++) msg[i] = '\0';
02746 
02747    if (!strcmp(param, "AZ_none")) {
02748       aztec_options_[AZ_scaling] = AZ_none;
02749       sprintf(msg, "No scaling");
02750    }
02751    else if (!strcmp(param, "AZ_Jacobi")) {
02752       aztec_options_[AZ_scaling] = AZ_Jacobi;
02753       sprintf(msg, "AZ_Jacobi scaling");
02754    }
02755    else if (!strcmp(param, "AZ_BJacobi")) {
02756       aztec_options_[AZ_scaling] = AZ_BJacobi;
02757       sprintf(msg, "AZ_BJacobi scaling");
02758    }
02759    else if (!strcmp(param, "AZ_row_sum")) {
02760       aztec_options_[AZ_scaling] = AZ_row_sum;
02761       sprintf(msg, "AZ_row_sum scaling");
02762    }
02763    else if (!strcmp(param, "AZ_sym_diag")) {
02764       aztec_options_[AZ_scaling] = AZ_sym_diag;
02765       sprintf(msg, "AZ_sym_diag scaling");
02766    }
02767    else if (!strcmp(param, "AZ_sym_row_sum")) {
02768       aztec_options_[AZ_scaling] = AZ_sym_row_sum;
02769       sprintf(msg, "AZ_sym_row_sum scaling");
02770    }
02771    else {
02772       if (thisProc_ == 0) {
02773          FEI_COUT << "Aztec_LinSysCore: Warning: requested scaling <" << param << "> not recognized." << FEI_ENDL;
02774       }
02775    }
02776 
02777    debugOutput(msg);
02778 
02779    delete [] msg;
02780 
02781    return;
02782 }
02783 
02784 //==============================================================================
02785 void Aztec_LinSysCore::setConvTest(const char* param) {
02786 
02787    char msg[64];
02788 
02789    if (!strcmp(param, "AZ_r0")) {
02790       aztec_options_[AZ_conv] = AZ_r0;
02791       sprintf(msg, "AZ_conv AZ_r0");
02792    }
02793    else if (!strcmp(param, "AZ_rhs")) {
02794       aztec_options_[AZ_conv] = AZ_rhs;
02795       sprintf(msg, "AZ_conv AZ_rhs");
02796    }
02797    else if (!strcmp(param, "AZ_Anorm")) {
02798       aztec_options_[AZ_conv] = AZ_Anorm;
02799       sprintf(msg, "AZ_conv AZ_Anorm");
02800    }
02801    else if (!strcmp(param, "AZ_sol")) {
02802       aztec_options_[AZ_conv] = AZ_sol;
02803       sprintf(msg, "AZ_conv AZ_sol");
02804    }
02805    else if (!strcmp(param, "AZ_weighted")) {
02806       aztec_options_[AZ_conv] = AZ_weighted;
02807       sprintf(msg, "AZ_conv AZ_weighted");
02808    }
02809    else if (!strcmp(param, "AZ_noscaled")) {
02810       aztec_options_[AZ_conv] = AZ_noscaled;
02811       sprintf(msg, "AZ_conv AZ_noscaled");
02812    }
02813    else {
02814       if (thisProc_ == 0) {
02815          FEI_COUT << "Aztec_LinSysCore: Warning: requested convergence test <" << param << "> not recognized." << FEI_ENDL;
02816       }
02817    }
02818 
02819    debugOutput(msg);
02820 
02821     return;
02822 }
02823 
02824 //==============================================================================
02825 void Aztec_LinSysCore::setPreCalc(const char* param)
02826 {
02827    char msg[64];
02828 
02829    if (!strcmp(param, "AZ_calc")) {
02830       aztec_options_[AZ_pre_calc] = AZ_calc;
02831       sprintf(msg, "AZ_pre_calc AZ_calc");
02832    }
02833    else if (!strcmp(param, "AZ_recalc")) {
02834       aztec_options_[AZ_pre_calc] = AZ_recalc;
02835       sprintf(msg, "AZ_pre_calc AZ_recalc");
02836    }
02837    else if (!strcmp(param, "AZ_reuse")) {
02838       aztec_options_[AZ_pre_calc] = AZ_reuse;
02839       sprintf(msg, "AZ_pre_calc AZ_reuse");
02840    }
02841    else {
02842       if (thisProc_ == 0) {
02843          FEI_COUT << "Aztec_LinSysCore: Warning: requested pre_calc <" << param << "> not recognized." << FEI_ENDL;
02844       }
02845    }
02846 
02847    debugOutput(msg);
02848 
02849    return;
02850 }
02851 
02852 //==============================================================================
02853 void Aztec_LinSysCore::setOverlap(const char* param)
02854 {
02855    char msg[64];
02856 
02857    if (!strcmp(param, "AZ_none")) {
02858       aztec_options_[AZ_overlap] = AZ_none;
02859       sprintf(msg, "AZ_overlap AZ_none");
02860    }
02861    else if (!strcmp(param, "AZ_diag")) {
02862       aztec_options_[AZ_overlap] = AZ_diag;
02863       sprintf(msg, "AZ_overlap AZ_diag");
02864    }
02865    else if (!strcmp(param, "AZ_full")) {
02866       aztec_options_[AZ_overlap] = AZ_full;
02867       sprintf(msg, "AZ_overlap AZ_full");
02868    }
02869    else {
02870       checkForOption("AZ_overlap", numParams_, paramStrings_,
02871                      aztec_options_[AZ_overlap]);
02872    }
02873 
02874    debugOutput(msg);
02875 
02876    return;
02877 }
02878 
02879 //==============================================================================
02880 void Aztec_LinSysCore::setOrthog(const char* param)
02881 {
02882    char msg[64];
02883 
02884    if (!strcmp(param, "AZ_classic")) {
02885       aztec_options_[AZ_orthog] = AZ_classic;
02886       sprintf(msg, "AZ_orthog AZ_classic");
02887    }
02888    else if (!strcmp(param, "AZ_modified")) {
02889       aztec_options_[AZ_orthog] = AZ_modified;
02890       sprintf(msg, "AZ_orthog AZ_modified");
02891    }
02892    else {
02893       if (thisProc_ == 0) {
02894          FEI_COUT << "Aztec_LinSysCore: Warning: requested orthog. <" << param << "> not recognized." << FEI_ENDL;
02895       }
02896    }
02897 
02898    debugOutput(msg);
02899 
02900    return;
02901 }
02902 
02903 //==============================================================================
02904 void Aztec_LinSysCore::setAuxVec(const char* param)
02905 {
02906    char msg[64];
02907 
02908    if (!strcmp(param, "AZ_resid")) {
02909       aztec_options_[AZ_aux_vec] = AZ_resid;
02910       sprintf(msg, "AZ_aux_vec AZ_resid");
02911    }
02912    else if (!strcmp(param, "AZ_rand")) {
02913       aztec_options_[AZ_aux_vec] = AZ_rand;
02914       sprintf(msg, "AZ_aux_vec AZ_rand");
02915    }
02916    else {
02917       if (thisProc_ == 0) {
02918          FEI_COUT << "Aztec_LinSysCore: Warning: requested aux_vec <" << param << "> not recognized." << FEI_ENDL;
02919       }
02920    }
02921 
02922    debugOutput(msg);
02923 
02924    return;
02925 }
02926 
02927 //==============================================================================
02928 void Aztec_LinSysCore::setAZ_output(const char* param)
02929 {
02930    char msg[64];
02931    int out = -1;
02932    int num = sscanf(param, "%d", &out);
02933    if (num == 1 && out > -1) {
02934      sprintf(msg, "AZ_output %d", out);
02935      aztec_options_[AZ_output] = out;
02936    }
02937    else if (!strcmp(param, "AZ_all")) {
02938       aztec_options_[AZ_output] = AZ_all;
02939       sprintf(msg, "AZ_output AZ_all");
02940    }
02941    else if (!strcmp(param, "AZ_none")) {
02942       aztec_options_[AZ_output] = AZ_none;
02943       sprintf(msg, "AZ_output AZ_none");
02944    }
02945    else if (!strcmp(param, "AZ_warnings")) {
02946       aztec_options_[AZ_output] = AZ_warnings;
02947       sprintf(msg, "AZ_output AZ_warnings");
02948    }
02949    else if (!strcmp(param, "AZ_last")) {
02950       aztec_options_[AZ_output] = AZ_last;
02951       sprintf(msg, "AZ_output AZ_last");
02952    }
02953    else {
02954       if (thisProc_ == 0) {
02955          FEI_COUT << "Aztec_LinSysCore: Warning: requested AZ_output <" << param << "> not recognized." << FEI_ENDL;
02956       }
02957    }
02958 
02959    debugOutput(msg);
02960 
02961    return;
02962 }
02963 
02964 //==============================================================================
02965 void Aztec_LinSysCore::checkForParam(const char* paramName,
02966                                      int numParams, char** paramStrings,
02967                                      double& param) {
02968    const char* parameter =
02969      snl_fei::getParamValue(paramName, numParams, paramStrings);
02970    if (parameter != NULL) {
02971       sscanf(parameter, "%le", &param);
02972    }
02973 }
02974 
02975 //==============================================================================
02976 void Aztec_LinSysCore::checkForOption(const char* paramName,
02977                                       int numParams, char** paramStrings,
02978                                       int& param) {
02979   const char* parameter =
02980     snl_fei::getParamValue(paramName, numParams, paramStrings);
02981    if (parameter != NULL) {
02982       sscanf(parameter, "%d", &param);
02983    }
02984 }
02985 
02986 //==============================================================================
02987 void Aztec_LinSysCore::setDebugOutput(const char* path, const char* name){
02988 //
02989 //This function turns on debug output, and opens a file to put it in.
02990 //
02991    if (debugOutput_) {
02992       fprintf(debugFile_,"setDebugOutput closing this file.");
02993       fflush(debugFile_);
02994       fclose(debugFile_);
02995       debugFile_ = NULL;
02996    }
02997 
02998    int pathLength = strlen(path);
02999    if (path != debugPath_) {
03000       delete [] debugPath_;
03001       debugPath_ = new char[pathLength + 1];
03002       sprintf(debugPath_, path);
03003    }
03004 
03005    int nameLength = strlen(name);
03006    if (name != debugFileName_) {
03007       delete [] debugFileName_;
03008       debugFileName_ = new char[nameLength + 1];
03009       sprintf(debugFileName_,name);
03010    }
03011 
03012    char* dbFileName = new char[pathLength + nameLength + 3];
03013 
03014    sprintf(dbFileName, "%s/%s", path, name);
03015 
03016    debugOutput_ = 1;
03017    debugFile_ = fopen(dbFileName, "w");
03018 
03019    if (!debugFile_){
03020       fei::console_out() << "couldn't open debug output file: " << dbFileName << FEI_ENDL;
03021       debugOutput_ = 0;
03022       delete [] debugPath_;
03023       debugPath_ = NULL;
03024       delete [] debugFileName_;
03025       debugFileName_ = NULL;
03026    }
03027 
03028    delete [] dbFileName;
03029 }
03030 
03031 //==============================================================================
03032 int Aztec_LinSysCore::VBRmatPlusScaledMat(AztecDVBR_Matrix* A,
03033                                            double scalar,
03034                                            AztecDVBR_Matrix* source)
03035 {
03036    int* nnz = new int[numLocalEqnBlks_];
03037    int* nblk = new int[numLocalEqnBlks_];
03038    int* src_nnz = new int[numLocalEqnBlks_];
03039    int* src_nblk = new int[numLocalEqnBlks_];
03040 
03041    if (nnz == NULL || nblk == NULL || src_nnz==NULL || src_nblk==NULL) {
03042       messageAbort("VBRMatPlusScaledMat: allocation failed");
03043    }
03044 
03045    A->getNumNonzerosPerRow(nnz);
03046    A->getNumBlocksPerRow(nblk);
03047    source->getNumNonzerosPerRow(src_nnz);
03048    source->getNumBlocksPerRow(src_nblk);
03049 
03050    int i, max_nnz = 0, max_nblk = 0;
03051    for(i=0; i<numLocalEqnBlks_; i++) {
03052       if (nnz[i] != src_nnz[i] || nblk[i] != src_nblk[i]) {
03053          messageAbort("VBRmatPlusScaledMat: matrix sizes don't match.");
03054       }
03055       if (max_nnz < nnz[i]) max_nnz = nnz[i];
03056       if (max_nblk < nblk[i]) max_nblk = nblk[i];
03057    }
03058    
03059    delete [] nnz;
03060    delete [] nblk;
03061    delete [] src_nnz;
03062    delete [] src_nblk;
03063 
03064    double* val = new double[max_nnz];
03065    int* colInds = new int[max_nblk];
03066    if (val==NULL || colInds==NULL) {
03067       messageAbort("VBRmatPlusScaledMat: allocation failed");
03068    }
03069    int len, nnzBlks;
03070 
03071    for(i=0; i<numLocalEqnBlks_; i++) {
03072       int row = localBlkOffset_+i;
03073       int err = source->getNumBlocksPerRow(row, nnzBlks);
03074       err += source->getNumNonzerosPerRow(row, len);
03075       err += source->getBlockRow(row, val, colInds, nnzBlks);
03076 
03077       if (err) messageAbort("VBRmatPlusScaledMat: error getting src row");
03078 
03079       for(int j=0; j<len; j++) val[j] *= scalar;
03080 
03081       err = A->sumIntoBlockRow(row, val, colInds, nnzBlks);
03082       if (err) messageAbort("VBRmatPlusScaledMat: error summing in row");
03083    }
03084 
03085    delete [] val;
03086    delete [] colInds;
03087    return(0);
03088 }
03089 
03090 //==============================================================================
03091 int Aztec_LinSysCore::MSRmatPlusScaledMat(AztecDMSR_Matrix* A,
03092                                            double scalar,
03093                                            AztecDMSR_Matrix* source)
03094 {
03095   return(A->addScaledMatrix(scalar, *source));
03096 }
03097 
03098 //==============================================================================
03099 void Aztec_LinSysCore::debugOutput(const char* msg) const {
03100    if (debugOutput_) {
03101       fprintf(debugFile_, "%s\n", msg);
03102       fflush(debugFile_);
03103    }
03104 }
03105 
03106 //==============================================================================
03107 int Aztec_LinSysCore::messageAbort(const char* msg) const {
03108    fei::console_out() << "Aztec_LinSysCore: " << msg << " Aborting." << FEI_ENDL;
03109 #ifndef FEI_SER
03110    MPI_Abort(comm_, -1);
03111 #else
03112    abort();
03113 #endif
03114    return(-1);
03115 }
03116 
03117 }//namespace fei_trilinos
03118 
03119 #endif
03120 //HAVE_FEI_AZTECOO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends