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