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