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

Generated on Wed May 12 21:30:40 2010 for FEI by  doxygen 1.4.7