FEI Version of the Day
snl_fei_LinearSystem_General.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_macros.hpp>
00010 #include <fei_utils.hpp>
00011 
00012 #include <cmath>
00013 
00014 #include <snl_fei_LinearSystem_General.hpp>
00015 #include <fei_MatrixReducer.hpp>
00016 #include <fei_Matrix_Impl.hpp>
00017 #include <fei_VectorSpace.hpp>
00018 #include <fei_MatrixGraph.hpp>
00019 #include <fei_SparseRowGraph.hpp>
00020 #include <snl_fei_Constraint.hpp>
00021 #include <fei_Record.hpp>
00022 #include <fei_utils.hpp>
00023 #include <fei_impl_utils.hpp>
00024 #include <fei_LogManager.hpp>
00025 
00026 #include <fei_DirichletBCRecord.hpp>
00027 #include <fei_DirichletBCManager.hpp>
00028 #include <fei_EqnBuffer.hpp>
00029 #include <fei_LinSysCoreFilter.hpp>
00030 
00031 #undef fei_file
00032 #define fei_file "snl_fei_LinearSystem_General.cpp"
00033 #include <fei_ErrMacros.hpp>
00034 
00035 //----------------------------------------------------------------------------
00036 snl_fei::LinearSystem_General::LinearSystem_General(fei::SharedPtr<fei::MatrixGraph>& matrixGraph)
00037   : fei::LinearSystem(matrixGraph),
00038     comm_(matrixGraph->getRowSpace()->getCommunicator()),
00039     essBCvalues_(NULL),
00040     resolveConflictRequested_(false),
00041     bcs_trump_slaves_(false),
00042     explicitBCenforcement_(false),
00043     BCenforcement_no_column_mod_(false),
00044     localProc_(0),
00045     numProcs_(1),
00046     name_(),
00047     named_loadcomplete_counter_(),
00048     iwork_(),
00049     dwork_(),
00050     dbgprefix_("LinSysG: ")
00051 {
00052   localProc_ = fei::localProc(comm_);
00053   numProcs_  = fei::numProcs(comm_);
00054 
00055   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
00056 
00057   std::vector<int> offsets;
00058   vecSpace->getGlobalIndexOffsets(offsets);
00059 
00060   firstLocalOffset_ = offsets[localProc_];
00061   lastLocalOffset_ = offsets[localProc_+1]-1;
00062 
00063   setName("dbg");
00064 }
00065 
00066 //----------------------------------------------------------------------------
00067 snl_fei::LinearSystem_General::~LinearSystem_General()
00068 {
00069   delete essBCvalues_;
00070 }
00071 
00072 //----------------------------------------------------------------------------
00073 int snl_fei::LinearSystem_General::parameters(int numParams,
00074            const char* const* paramStrings)
00075 {
00076   if (numParams == 0 || paramStrings == NULL) return(0);
00077 
00078   const char* param = snl_fei::getParam("name", numParams, paramStrings);
00079   if (param != NULL) {
00080     if (strlen(param) < 6) ERReturn(-1);
00081 
00082     setName(&(param[5]));
00083   }
00084 
00085   param = snl_fei::getParam("resolveConflict",numParams,paramStrings);
00086   if (param != NULL){
00087     resolveConflictRequested_ = true;
00088   }
00089 
00090   param = snl_fei::getParam("BCS_TRUMP_SLAVE_CONSTRAINTS",
00091                             numParams,paramStrings);
00092   if (param != NULL) {
00093     bcs_trump_slaves_ = true;
00094   }
00095 
00096   param = snl_fei::getParam("EXPLICIT_BC_ENFORCEMENT",numParams,paramStrings);
00097   if (param != NULL){
00098     explicitBCenforcement_ = true;
00099   }
00100 
00101   param = snl_fei::getParam("BC_ENFORCEMENT_NO_COLUMN_MOD",numParams,paramStrings);
00102   if (param != NULL){
00103     BCenforcement_no_column_mod_ = true;
00104   }
00105 
00106   param = snl_fei::getParamValue("FEI_OUTPUT_LEVEL",numParams,paramStrings);
00107   if (param != NULL) {
00108     setOutputLevel(fei::utils::string_to_output_level(param));
00109   }
00110 
00111   if (matrix_.get() != NULL) {
00112     fei::Matrix* matptr = matrix_.get();
00113     fei::MatrixReducer* matred = dynamic_cast<fei::MatrixReducer*>(matptr);
00114     if (matred != NULL) {
00115       matptr = matred->getTargetMatrix().get();
00116     }
00117     fei::Matrix_Impl<LinearSystemCore>* lscmatrix =
00118       dynamic_cast<fei::Matrix_Impl<LinearSystemCore>*>(matptr);
00119     if (lscmatrix != NULL) {
00120       lscmatrix->getMatrix()->parameters(numParams, (char**)paramStrings);
00121     }
00122   }
00123 
00124   return(0);
00125 }
00126 
00127 //----------------------------------------------------------------------------
00128 int snl_fei::LinearSystem_General::parameters(const fei::ParameterSet& params)
00129 {
00130   int numParams = 0;
00131   const char** paramStrings = NULL;
00132   std::vector<std::string> stdstrings;
00133   fei::utils::convert_ParameterSet_to_strings(&params, stdstrings);
00134   fei::utils::strings_to_char_ptrs(stdstrings, numParams, paramStrings);
00135 
00136   int err = parameters(numParams, paramStrings);
00137 
00138   delete [] paramStrings;
00139 
00140   return(err);
00141 }
00142 
00143 //----------------------------------------------------------------------------
00144 void snl_fei::LinearSystem_General::setName(const char* name)
00145 {
00146   if (name == NULL) return;
00147 
00148   if (name_ == name) return;
00149 
00150   name_ = name;
00151 
00152   std::map<std::string,unsigned>::iterator
00153     iter = named_loadcomplete_counter_.find(name_);
00154   if (iter == named_loadcomplete_counter_.end()) {
00155     static int counter = 0;
00156     named_loadcomplete_counter_.insert(std::make_pair(name_, counter));
00157     ++counter;
00158   }
00159 }
00160 
00161 //----------------------------------------------------------------------------
00162 int snl_fei::LinearSystem_General::loadEssentialBCs(int numIDs,
00163                          const int* IDs,
00164                          int idType,
00165                          int fieldID,
00166                          int offsetIntoField,
00167                          const double* prescribedValues)
00168 {
00169   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != 0) {
00170     FEI_OSTREAM& os = *output_stream_;
00171     os << "loadEssentialBCs, numIDs: "<<numIDs<<", idType: " <<idType
00172     <<", fieldID: "<<fieldID<<", offsetIntoField: "<<offsetIntoField<<FEI_ENDL;
00173   }
00174 
00175   return fei::LinearSystem::loadEssentialBCs(numIDs, IDs, idType, fieldID,
00176                                            offsetIntoField, prescribedValues);
00177 }
00178 
00179 //----------------------------------------------------------------------------
00180 int snl_fei::LinearSystem_General::loadEssentialBCs(int numIDs,
00181                          const int* IDs,
00182                          int idType,
00183                          int fieldID,
00184                          const int* offsetsIntoField,
00185                          const double* prescribedValues)
00186 {
00187   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != 0) {
00188     FEI_OSTREAM& os = *output_stream_;
00189     for(int i=0; i<numIDs; ++i) {
00190       os << "loadEssentialBCs idType: " <<idType
00191         <<", fieldID: "<<fieldID<<", ID: " << IDs[i]<<", offsetIntoField: "<<offsetsIntoField[i]<<", val: " << prescribedValues[i] << FEI_ENDL;
00192     }
00193   }
00194 
00195   return fei::LinearSystem::loadEssentialBCs(numIDs, IDs, idType, fieldID,
00196                                            offsetsIntoField, prescribedValues);
00197 }
00198 
00199 //----------------------------------------------------------------------------
00200 int snl_fei::LinearSystem_General::loadComplete(bool applyBCs,
00201                                                 bool globalAssemble)
00202 {
00203   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != 0) {
00204     FEI_OSTREAM& os = *output_stream_;
00205     os << dbgprefix_<<"loadComplete"<<FEI_ENDL;
00206   }
00207 
00208   if (dbcManager_ == NULL) {
00209     dbcManager_ = new fei::DirichletBCManager(matrixGraph_->getRowSpace());
00210   }
00211 
00212   if (globalAssemble) {
00213 
00214     if (matrix_.get() != NULL) {
00215       CHK_ERR( matrix_->gatherFromOverlap() );
00216     }
00217 
00218     if (rhs_.get() != NULL) {
00219       CHK_ERR( rhs_->gatherFromOverlap() );
00220     }
00221 
00222   }
00223 
00224   unsigned counter = 0;
00225 
00226   std::map<std::string,unsigned>::iterator
00227     iter = named_loadcomplete_counter_.find(name_);
00228   if (iter == named_loadcomplete_counter_.end()) {
00229     FEI_COUT << "fei::LinearSystem::loadComplete internal error, name "
00230       << name_ << " not found." << FEI_ENDL;
00231   }
00232   else {
00233     counter = iter->second++;
00234   }
00235 
00236   if (output_level_ >= fei::FULL_LOGS) {
00237     std::string opath = fei::LogManager::getLogManager().getOutputPath();
00238     if (opath == "") opath = ".";
00239 
00240     FEI_OSTRINGSTREAM Aname;
00241     FEI_OSTRINGSTREAM bname;
00242     FEI_OSTRINGSTREAM xname;
00243     Aname << opath << "/";
00244     bname << opath << "/";
00245     xname << opath << "/";
00246 
00247     Aname << "A_"<<name_<<".preBC.np"<<numProcs_<<".slv"<<counter<< ".mtx";
00248 
00249     bname << "b_"<<name_<<".preBC.np"<<numProcs_<<".slv"<<counter<< ".vec";
00250 
00251     std::string Aname_str = Aname.str();
00252     const char* Aname_c_str = Aname_str.c_str();
00253     CHK_ERR( matrix_->writeToFile(Aname_c_str) );
00254 
00255     std::string bname_str = bname.str();
00256     const char* bname_c_str = bname_str.c_str();
00257     CHK_ERR( rhs_->writeToFile(bname_c_str) );
00258   }
00259 
00260   CHK_ERR( implementBCs(applyBCs) );
00261 
00262   if (globalAssemble) {
00263     CHK_ERR( matrix_->globalAssemble() );
00264   }
00265 
00266   if (output_level_ == fei::STATS || output_level_ == fei::ALL) {
00267     int globalNumSlaveCRs = matrixGraph_->getGlobalNumSlaveConstraints();
00268     if (localProc_ == 0) {
00269       FEI_COUT << "Global Neqns: " << matrix_->getGlobalNumRows();
00270       if (globalNumSlaveCRs > 0) {
00271   FEI_COUT << ", Global NslaveCRs: " << globalNumSlaveCRs;
00272       }
00273       FEI_COUT << FEI_ENDL;
00274     }
00275   }
00276 
00277   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00278     FEI_OSTREAM& os = *output_stream_;
00279     os << dbgprefix_<<"Neqns=" << matrix_->getGlobalNumRows();
00280     int globalNumSlaveCRs = matrixGraph_->getGlobalNumSlaveConstraints();
00281     if (globalNumSlaveCRs > 0) {
00282       os << ", Global NslaveCRs=" << globalNumSlaveCRs;
00283     }
00284     os << FEI_ENDL;
00285   }
00286 
00287   if (output_level_ >= fei::MATRIX_FILES) {
00288     std::string opath = fei::LogManager::getLogManager().getOutputPath();
00289     if (opath == "") opath = ".";
00290 
00291     FEI_OSTRINGSTREAM Aname;
00292     FEI_OSTRINGSTREAM bname;
00293     FEI_OSTRINGSTREAM xname;
00294     Aname << opath << "/";
00295     bname << opath << "/";
00296     xname << opath << "/";
00297 
00298     Aname << "A_" <<name_<<".np"<<numProcs_<< ".slv" << counter << ".mtx";
00299 
00300     bname << "b_" <<name_<<".np"<<numProcs_<< ".slv" << counter << ".vec";
00301 
00302     xname << "x0_" <<name_<<".np"<<numProcs_<< ".slv" << counter << ".vec";
00303 
00304     std::string Aname_str = Aname.str();
00305     const char* Aname_c_str = Aname_str.c_str();
00306     CHK_ERR( matrix_->writeToFile(Aname_c_str) );
00307 
00308     std::string bname_str = bname.str();
00309     const char* bname_c_str = bname_str.c_str();
00310     CHK_ERR( rhs_->writeToFile(bname_c_str) );
00311 
00312     std::string xname_str = xname.str();
00313     const char* xname_c_str = xname_str.c_str();
00314     CHK_ERR( soln_->writeToFile(xname_c_str) );
00315   }
00316 
00317   return(0);
00318 }
00319 
00320 //----------------------------------------------------------------------------
00321 int snl_fei::LinearSystem_General::setBCValuesOnVector(fei::Vector* vector)
00322 {
00323   if (essBCvalues_ == NULL) {
00324     return(0);
00325   }
00326 
00327   if (essBCvalues_->size() == 0) {
00328     return(0);
00329   }
00330 
00331   CHK_ERR( vector->copyIn(essBCvalues_->size(),
00332         &(essBCvalues_->indices())[0],
00333         &(essBCvalues_->coefs())[0]) );
00334 
00335   return(0);
00336 }
00337 
00338 //----------------------------------------------------------------------------
00339 bool snl_fei::LinearSystem_General::eqnIsEssentialBC(int globalEqnIndex) const
00340 {
00341   if (essBCvalues_ == NULL) return(false);
00342 
00343   std::vector<int>& indices = essBCvalues_->indices();
00344   int offset = fei::binarySearch(globalEqnIndex, indices);
00345   return( offset < 0 ? false : true);
00346 }
00347 
00348 //----------------------------------------------------------------------------
00349 void snl_fei::LinearSystem_General::getEssentialBCs(std::vector<int>& bcEqns,
00350                                              std::vector<double>& bcVals) const
00351 {
00352   bcEqns.clear();
00353   bcVals.clear();
00354   if (essBCvalues_ == NULL) return;
00355 
00356   int num = essBCvalues_->size();
00357   bcEqns.resize(num);
00358   bcVals.resize(num);
00359   int* essbcs = &(essBCvalues_->indices())[0];
00360   double* vals = &(essBCvalues_->coefs())[0];
00361   for(int i=0; i<num; ++i) {
00362     bcEqns[i] = essbcs[i];
00363     bcVals[i] = vals[i];
00364   }
00365 }
00366 
00367 //----------------------------------------------------------------------------
00368 void snl_fei::LinearSystem_General::getConstrainedEqns(std::vector<int>& crEqns) const
00369 {
00370   matrixGraph_->getConstrainedIndices(crEqns);
00371 }
00372 
00373 //----------------------------------------------------------------------------
00374 int extractDirichletBCs(fei::DirichletBCManager* bcManager,
00375                 fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00376                 fei::CSVec* essBCvalues,
00377                 bool resolveConflictRequested,
00378                 bool bcs_trump_slaves)
00379 {
00380 //  int numLocalBCs = bcManager->getNumBCRecords();
00381 //  int globalNumBCs = 0;
00382 //  MPI_Comm comm = matrixGraph->getRowSpace()->getCommunicator();
00383 //  fei::GlobalSum(comm, numLocalBCs, globalNumBCs);
00384 //  if (globalNumBCs == 0) {
00385 //    return(0);
00386 //  }
00387 
00388   fei::SharedPtr<fei::FillableMat> localBCeqns(new fei::FillableMat);
00389   fei::SharedPtr<fei::Matrix_Impl<fei::FillableMat> > bcEqns;
00390 //  matrixGraph->getRowSpace()->initComplete();
00391   int numSlaves = matrixGraph->getGlobalNumSlaveConstraints();
00392   fei::SharedPtr<fei::Reducer> reducer = matrixGraph->getReducer();
00393 
00394   int numIndices = numSlaves>0 ?
00395     reducer->getLocalReducedEqns().size() :
00396     matrixGraph->getRowSpace()->getNumIndices_Owned();
00397 
00398   bool zeroSharedRows = false;
00399   bcEqns.reset(new fei::Matrix_Impl<fei::FillableMat>(localBCeqns, matrixGraph, numIndices, zeroSharedRows));
00400   fei::SharedPtr<fei::Matrix> bcEqns_reducer;
00401   if (numSlaves > 0) {
00402     bcEqns_reducer.reset(new fei::MatrixReducer(reducer, bcEqns));
00403   }
00404 
00405   fei::Matrix& bcEqns_mat = bcEqns_reducer.get()==NULL ?
00406       *bcEqns : *bcEqns_reducer;
00407 
00408   CHK_ERR( bcManager->finalizeBCEqns(bcEqns_mat, bcs_trump_slaves) );
00409 
00410   if (resolveConflictRequested) {
00411     fei::SharedPtr<fei::FillableMat> mat = bcEqns->getMatrix();
00412     std::vector<int> bcEqnNumbers;
00413     fei::get_row_numbers(*mat, bcEqnNumbers);
00414     CHK_ERR( snl_fei::resolveConflictingCRs(*matrixGraph, bcEqns_mat,
00415                                             bcEqnNumbers) );
00416   }
00417 
00418   std::vector<int> essEqns;
00419   std::vector<double> values;
00420 
00421   std::map<int,fei::FillableMat*>& remotes = bcEqns->getRemotelyOwnedMatrices();
00422   std::map<int,fei::FillableMat*>::iterator
00423     it = remotes.begin(),
00424     it_end = remotes.end();
00425   for(; it!=it_end; ++it) {
00426     fei::impl_utils::separate_BC_eqns( *(it->second), essEqns, values);
00427   }
00428 
00429 //  CHK_ERR( bcEqns->gatherFromOverlap(false) );
00430 
00431   fei::impl_utils::separate_BC_eqns( *(bcEqns->getMatrix()), essEqns, values);
00432 
00433   if (essEqns.size() > 0) {
00434     int* essEqnsPtr = &essEqns[0];
00435     double* valuesPtr = &values[0];
00436 
00437     for(unsigned i=0; i<essEqns.size(); ++i) {
00438       int eqn = essEqnsPtr[i];
00439       double value = valuesPtr[i];
00440       fei::put_entry(*essBCvalues, eqn, value);
00441     }
00442   }
00443 
00444   return(0);
00445 }
00446 
00447 //----------------------------------------------------------------------------
00448 int snl_fei::LinearSystem_General::implementBCs(bool applyBCs)
00449 {
00450   if (essBCvalues_ != NULL) {
00451     delete essBCvalues_;
00452   }
00453 
00454   essBCvalues_ = new fei::CSVec;
00455 
00456   CHK_ERR( extractDirichletBCs(dbcManager_, matrixGraph_,
00457                        essBCvalues_,  resolveConflictRequested_,
00458                       bcs_trump_slaves_) );
00459 
00460   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00461     FEI_OSTREAM& os = *output_stream_;
00462     std::vector<int>& indices = essBCvalues_->indices();
00463     std::vector<double>& coefs= essBCvalues_->coefs();
00464     for(size_t i=0; i<essBCvalues_->size(); ++i) {
00465       os << "essBCeqns["<<i<<"]: "<<indices[i]<<", "<<coefs[i]<<FEI_ENDL;
00466     }
00467   }
00468 
00469   //If the underlying matrix is a LinearSystemCore instance, then this
00470   //function will return 0, and we're done. A non-zero return-code means
00471   //we should continue and enforce the BCs assuming a general matrix.
00472 
00473   int returncode = enforceEssentialBC_LinSysCore();
00474   if (returncode == 0) {
00475     return(0);
00476   }
00477 
00478   fei::CSVec allEssBCs;
00479   if (!BCenforcement_no_column_mod_) {
00480     fei::impl_utils::global_union(comm_, *essBCvalues_, allEssBCs);
00481 
00482     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00483       FEI_OSTREAM& os = *output_stream_;
00484       os << "  implementBCs, essBCvalues_.size(): "<<essBCvalues_->size()
00485          << ", allEssBCs.size(): " << allEssBCs.size()<<FEI_ENDL;
00486     }
00487   }
00488 
00489   if (essBCvalues_->size() > 0) {
00490     enforceEssentialBC_step_1(*essBCvalues_);
00491   }
00492 
00493   if (!BCenforcement_no_column_mod_ && allEssBCs.size() > 0) {
00494     enforceEssentialBC_step_2(allEssBCs);
00495   }
00496 
00497   return(0);
00498 }
00499 
00500 //----------------------------------------------------------------------------
00501 int snl_fei::LinearSystem_General::enforceEssentialBC_LinSysCore()
00502 {
00503   fei::Matrix* matptr = matrix_.get();
00504   fei::MatrixReducer* matred = dynamic_cast<fei::MatrixReducer*>(matptr);
00505   if (matred != NULL) {
00506     matptr = matred->getTargetMatrix().get();
00507   }
00508 
00509   fei::Matrix_Impl<LinearSystemCore>* lscmatrix =
00510     dynamic_cast<fei::Matrix_Impl<LinearSystemCore>*>(matptr);
00511   if (lscmatrix == 0) {
00512     return(-1);
00513   }
00514 
00515   int localsize = matrixGraph_->getRowSpace()->getNumIndices_Owned();
00516   fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
00517   if (matrixGraph_->getGlobalNumSlaveConstraints() > 0) {
00518     localsize = reducer->getLocalReducedEqns().size();
00519   }
00520 
00521   fei::SharedPtr<fei::FillableMat> inner(new fei::FillableMat);
00522   bool zeroSharedRows = false;
00523   fei::SharedPtr<fei::Matrix_Impl<fei::FillableMat> > matrix;
00524   matrix.reset(new fei::Matrix_Impl<fei::FillableMat>(inner, matrixGraph_, localsize, zeroSharedRows));
00525 
00526   fei::SharedPtr<fei::SparseRowGraph> remoteGraph =
00527     matrixGraph_->getRemotelyOwnedGraphRows();
00528 
00529   if (!BCenforcement_no_column_mod_) {
00530     CHK_ERR( snl_fei::gatherRemoteEssBCs(*essBCvalues_, remoteGraph.get(), *matrix) );
00531   }
00532 
00533   unsigned numBCRows = inner->getNumRows();
00534 
00535   if (output_stream_ != NULL && output_level_ >= fei::BRIEF_LOGS) {
00536     FEI_OSTREAM& os = *output_stream_;
00537     os << "#enforceEssentialBC_LinSysCore RemEssBCs to enforce: "
00538        << numBCRows << FEI_ENDL;
00539   }
00540 
00541   if (numBCRows > 0 && !BCenforcement_no_column_mod_) {
00542     std::vector<int*> colIndices(numBCRows);
00543     std::vector<double*> coefs(numBCRows);
00544     std::vector<int> colIndLengths(numBCRows);
00545 
00546     fei::CSRMat csrmat(*inner);
00547     fei::SparseRowGraph& srg = csrmat.getGraph();
00548 
00549     int numEqns = csrmat.getNumRows();
00550     int* eqns = &(srg.rowNumbers[0]);
00551     int* rowOffsets = &(srg.rowOffsets[0]);
00552 
00553     for(int i=0; i<numEqns; ++i) {
00554       colIndices[i] = &(srg.packedColumnIndices[rowOffsets[i]]);
00555       coefs[i] = &(csrmat.getPackedCoefs()[rowOffsets[i]]);
00556       colIndLengths[i] = rowOffsets[i+1] - rowOffsets[i];
00557     }
00558 
00559     int** colInds = &colIndices[0];
00560     int* colIndLens = &colIndLengths[0];
00561     double** BCcoefs = &coefs[0];
00562 
00563     if (output_stream_ != NULL && output_level_ > fei::BRIEF_LOGS) {
00564       FEI_OSTREAM& os = *output_stream_;
00565       for(int i=0; i<numEqns; ++i) {
00566         os << "remBCeqn: " << eqns[i] << ", inds/coefs: ";
00567         for(int j=0; j<colIndLens[i]; ++j) {
00568           os << "("<<colInds[i][j]<<","<<BCcoefs[i][j]<<") ";
00569         }
00570         os << FEI_ENDL;
00571       }
00572     }
00573 
00574     int errcode = lscmatrix->getMatrix()->enforceRemoteEssBCs(numEqns,
00575                     eqns,
00576                     colInds,
00577                     colIndLens,
00578                     BCcoefs);
00579     if (errcode != 0) {
00580       return(errcode);
00581     }
00582   }
00583 
00584   int numEqns = essBCvalues_->size();
00585   if (numEqns > 0) {
00586     int* eqns = &(essBCvalues_->indices())[0];
00587     double* bccoefs = &(essBCvalues_->coefs())[0];
00588     std::vector<double> ones(numEqns, 1.0);
00589 
00590     return(lscmatrix->getMatrix()->enforceEssentialBC(eqns, &ones[0],
00591                 bccoefs, numEqns));
00592   }
00593 
00594   return(0);
00595 }
00596 
00597 //----------------------------------------------------------------------------
00598 void snl_fei::LinearSystem_General::enforceEssentialBC_step_1(fei::CSVec& essBCs)
00599 {
00600   //to enforce essential boundary conditions, we do the following:
00601   //
00602   //  1.  for each eqn (== essBCs.indices()[n]), {
00603   //        put zeros in row A[eqn], but leave 1.0 on the diagonal
00604   //        set b[eqn] = essBCs.coefs()[n]
00605   //      }
00606   //
00607   //  2.  for i in 1..numRows (i.e., all rows) {
00608   //        if (i in bcEqns) continue;
00609   //        b[i] -= A[i,eqn] * essBCs.coefs()[n]
00610   //        A[i,eqn] = 0.0;
00611   //      }
00612   //
00613   //It is important to note that for step 1, essBCs need only contain
00614   //local eqns, but for step 2 it should contain *ALL* bc eqns.
00615   //
00616   //This function performs step 1.
00617 
00618   int numEqns = essBCs.size();
00619   int* eqns = &(essBCs.indices())[0];
00620   double* bcCoefs = &(essBCs.coefs())[0];
00621 
00622   std::vector<double> coefs;
00623   std::vector<int> indices;
00624 
00625   fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
00626   bool haveSlaves = reducer.get()!=NULL;
00627 
00628   try {
00629   for(int i=0; i<numEqns; i++) {
00630     int eqn = eqns[i];
00631 
00632     //if slave-constraints are present, the incoming bc-eqns are in
00633     //the reduced equation space. so we actually have to translate them back
00634     //to the unreduced space before passing them into the fei::Matrix object,
00635     //because the fei::Matrix object has machinery to translate unreduced eqns
00636     //to the reduced space.
00637     //Also, our firstLocalOffset_ and lastLocalOffset_ attributes are in the
00638     //unreduced space.
00639     if (haveSlaves) {
00640       eqn = reducer->translateFromReducedEqn(eqn);
00641     }
00642 
00643     if (eqn < firstLocalOffset_ || eqn > lastLocalOffset_) continue;
00644 
00645     //put gamma/alpha on the rhs for this ess-BC equation.
00646     double bcValue = bcCoefs[i];
00647     int err = rhs_->copyIn(1, &eqn, &bcValue);
00648     if (err != 0) {
00649       FEI_OSTRINGSTREAM osstr;
00650       osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_1 ERROR: "
00651       << "err="<<err<<" returned from rhs_->copyIn row="<<eqn;
00652       throw std::runtime_error(osstr.str());
00653     }
00654 
00655     err = getMatrixRow(matrix_.get(), eqn, coefs, indices);
00656     if (err != 0 || indices.size() < 1) {
00657       continue;
00658     }
00659 
00660     int rowLen = indices.size();
00661     int* indPtr = &indices[0];
00662 
00663     //first, put zeros in the row and 1.0 on the diagonal...
00664     for(int j=0; j<rowLen; j++) {
00665       if (indPtr[j] == eqn) coefs[j] = 1.0;
00666       else coefs[j] = 0.0;
00667     }
00668 
00669     double* coefPtr = &coefs[0];
00670 
00671     err = matrix_->copyIn(1, &eqn, rowLen, indPtr, &coefPtr);
00672     if (err != 0) {
00673       FEI_OSTRINGSTREAM osstr;
00674       osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_1 ERROR: "
00675       << "err="<<err<<" returned from matrix_->copyIn row="<<eqn;
00676       throw std::runtime_error(osstr.str());
00677     }
00678   }//for i
00679   }
00680   catch(std::runtime_error& exc) {
00681     FEI_OSTRINGSTREAM osstr;
00682     osstr << "fei::LinearSystem::enforceEssentialBC: ERROR, caught exception: "
00683         << exc.what();
00684     throw std::runtime_error(osstr.str());
00685   }
00686 }
00687 
00688 //----------------------------------------------------------------------------
00689 void snl_fei::LinearSystem_General::enforceEssentialBC_step_2(fei::CSVec& essBCs)
00690 {
00691   //to enforce essential boundary conditions, we do the following:
00692   //
00693   //  1.  for each eqn (== essBCs.indices()[n]), {
00694   //        put zeros in row A[eqn], but leave 1.0 on the diagonal
00695   //        set b[eqn] = essBCs.coefs()[n]
00696   //      }
00697   //
00698   //  2.  for i in 1..numRows (i.e., all rows) {
00699   //        if (i in bcEqns) continue;
00700   //        b[i] -= A[i,eqn] * essBCs.coefs()[n]
00701   //        A[i,eqn] = 0.0;
00702   //      }
00703   //
00704   //It is important to note that for step 1, essBCs need only contain
00705   //local eqns, but for step 2 it should contain *ALL* bc eqns.
00706   //
00707   //This function performs step 2.
00708 
00709   int numBCeqns = essBCs.size();
00710   if (numBCeqns < 1) {
00711     return;
00712   }
00713 
00714   int* bcEqns = &(essBCs.indices())[0];
00715   double* bcCoefs = &(essBCs.coefs())[0];
00716 
00717   fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
00718   bool haveSlaves = reducer.get()!=NULL;
00719   if (haveSlaves) {
00720     for(int i=0; i<numBCeqns; ++i) {
00721       bcEqns[i] = reducer->translateFromReducedEqn(bcEqns[i]);
00722     }
00723   }
00724 
00725   int firstBCeqn = bcEqns[0];
00726   int lastBCeqn = bcEqns[numBCeqns-1];
00727 
00728   std::vector<double> coefs;
00729   std::vector<int> indices;
00730 
00731   int insertPoint;
00732 
00733   int nextBCeqnOffset = 0;
00734   int nextBCeqn = bcEqns[nextBCeqnOffset];
00735 
00736   for(int i=firstLocalOffset_; i<=lastLocalOffset_; ++i) {
00737     if (haveSlaves) {
00738       if (reducer->isSlaveEqn(i)) continue;
00739     }
00740 
00741     bool should_continue = false;
00742     if (i >= nextBCeqn) {
00743       if (i == nextBCeqn) {
00744   ++nextBCeqnOffset;
00745   if (nextBCeqnOffset < numBCeqns) {
00746     nextBCeqn = bcEqns[nextBCeqnOffset];
00747   }
00748   else {
00749     nextBCeqn = lastLocalOffset_+1;
00750   }
00751 
00752   should_continue = true;
00753       }
00754       else {
00755   while(nextBCeqn <= i) {
00756     if (nextBCeqn == i) should_continue = true;
00757     ++nextBCeqnOffset;
00758     if (nextBCeqnOffset < numBCeqns) {
00759       nextBCeqn = bcEqns[nextBCeqnOffset];
00760     }
00761     else {
00762       nextBCeqn = lastLocalOffset_+1;
00763     }
00764   }
00765       }
00766     }
00767 
00768     if (should_continue) continue;
00769 
00770     int err = getMatrixRow(matrix_.get(), i, coefs, indices);
00771     if (err != 0 || indices.size() < 1) {
00772       continue;
00773     }
00774 
00775     int numIndices = indices.size();
00776     int* indicesPtr = &indices[0];
00777     double* coefsPtr = &coefs[0];
00778     bool modifiedCoef = false;
00779 
00780     fei::insertion_sort_with_companions(numIndices, indicesPtr, coefsPtr);
00781 
00782     if (indicesPtr[0] > lastBCeqn || indicesPtr[numIndices-1] < firstBCeqn) {
00783       continue;
00784     }
00785 
00786     double value = 0.0;
00787     int offset = 0;
00788 
00789     for(int j=0; j<numIndices; ++j) {
00790       int idx = indicesPtr[j];
00791       offset = fei::binarySearch(idx, bcEqns, numBCeqns,
00792              insertPoint);
00793       if (offset > -1) {
00794   value -= bcCoefs[offset]*coefsPtr[j];
00795 
00796   coefsPtr[j] = 0.0;
00797   modifiedCoef = true;
00798       }
00799     }
00800 
00801     if (modifiedCoef) {
00802       err = matrix_->copyIn(1, &i, numIndices, indicesPtr, &coefsPtr);
00803       if (err != 0) {
00804   FEI_OSTRINGSTREAM osstr;
00805   osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_2 ERROR: "
00806         << "err="<<err<<" returned from matrix_->copyIn, row="<<i;
00807   throw std::runtime_error(osstr.str());
00808       }
00809     }
00810 
00811     const double fei_eps = 1.e-49;
00812     if (std::abs(value) > fei_eps) {
00813       rhs_->sumIn(1, &i, &value);
00814 
00815       if (output_level_ >= fei::FULL_LOGS && output_stream_ != 0) {
00816   FEI_OSTREAM& os = *output_stream_;
00817   os << "enfEssBC_step2: rhs["<<i<<"] += "<<value<<FEI_ENDL;
00818       }
00819     }
00820   }
00821 }
00822 
00823 //----------------------------------------------------------------------------
00824 int snl_fei::LinearSystem_General::getMatrixRow(fei::Matrix* matrix, int row,
00825             std::vector<double>& coefs,
00826             std::vector<int>& indices)
00827 {
00828   int len = 0;
00829   int err = matrix->getRowLength(row, len);
00830   if (err != 0) {
00831     coefs.resize(0);
00832     indices.resize(0);
00833     return(err);
00834   }
00835 
00836   if ((int)coefs.size() != len) {
00837     coefs.resize(len);
00838   }
00839   if ((int)indices.size() != len) {
00840     indices.resize(len);
00841   }
00842 
00843   CHK_ERR( matrix->copyOutRow(row, len, &coefs[0], &indices[0]));
00844 
00845   return(0);
00846 }
00847 
00848 //----------------------------------------------------------------------------
00849 int snl_fei::LinearSystem_General::loadLagrangeConstraint(int constraintID,
00850                 const double *weights,
00851                 double rhsValue)
00852 {
00853   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00854     FEI_OSTREAM& os = *output_stream_;
00855     os << "loadLagrangeConstraint crID: "<<constraintID<<FEI_ENDL;
00856   }
00857 
00858   Constraint<fei::Record<int>*>* cr =
00859     matrixGraph_->getLagrangeConstraint(constraintID);
00860   if (cr == NULL) {
00861     return(-1);
00862   }
00863 
00864   CHK_ERR( matrixGraph_->getConstraintConnectivityIndices(cr, iwork_) );
00865 
00866   //Let's attach the weights to the constraint-record now.
00867   std::vector<double>& cr_weights = cr->getMasterWeights();
00868   cr_weights.resize(iwork_.size());
00869   for(unsigned i=0; i<iwork_.size(); ++i) {
00870     cr_weights.push_back(weights[i]);
00871   }
00872 
00873   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph_->getRowSpace();
00874 
00875   int crEqn = -1;
00876   CHK_ERR( vecSpace->getGlobalIndex(cr->getIDType(),
00877              cr->getConstraintID(),
00878              crEqn) );
00879 
00880   //now add the row contribution to the matrix and rhs
00881   int numIndices = iwork_.size();
00882   int* indicesPtr = &(iwork_[0]);
00883 
00884   CHK_ERR( matrix_->sumIn(1, &crEqn, numIndices, indicesPtr, &weights) );
00885 
00886   CHK_ERR( rhs_->sumIn(1, &crEqn, &rhsValue) );
00887 
00888   //now add the column contributions to the matrix
00889   for(int k=0; k<numIndices; ++k) {
00890     double* thisWeight = (double*)(&(weights[k]));
00891     CHK_ERR( matrix_->sumIn(1, &(indicesPtr[k]), 1, &crEqn, &thisWeight) );
00892   }
00893 
00894   return(0);
00895 }
00896 
00897 //----------------------------------------------------------------------------
00898 int snl_fei::LinearSystem_General::loadPenaltyConstraint(int constraintID,
00899                const double *weights,
00900                double penaltyValue,
00901                double rhsValue)
00902 {
00903   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00904     FEI_OSTREAM& os = *output_stream_;
00905     os << "loadPenaltyConstraint crID: "<<constraintID<<FEI_ENDL;
00906   }
00907 
00908   Constraint<fei::Record<int>*>* cr =
00909     matrixGraph_->getPenaltyConstraint(constraintID);
00910   if (cr == NULL) {
00911     return(-1);
00912   }
00913 
00914   CHK_ERR( matrixGraph_->getConstraintConnectivityIndices(cr, iwork_) );
00915 
00916   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph_->getRowSpace();
00917 
00918   int numIndices = iwork_.size();
00919   int* indicesPtr = &(iwork_[0]);
00920 
00921   //now add the contributions to the matrix and rhs
00922   std::vector<double> coefs(numIndices);
00923   double* coefPtr = &coefs[0];
00924   for(int i=0; i<numIndices; ++i) {
00925     for(int j=0; j<numIndices; ++j) {
00926       coefPtr[j] = weights[i]*weights[j]*penaltyValue;
00927     }
00928     CHK_ERR( matrix_->sumIn(1, &(indicesPtr[i]), numIndices, indicesPtr,
00929           &coefPtr) );
00930 
00931     double rhsCoef = weights[i]*penaltyValue*rhsValue;
00932     CHK_ERR( rhs_->sumIn(1, &(indicesPtr[i]), &rhsCoef) );
00933   }
00934 
00935   return(0);
00936 }
00937 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends