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   bcEqns.reset(new fei::Matrix_Impl<fei::FillableMat>(localBCeqns, matrixGraph, numIndices));
00399   fei::SharedPtr<fei::Matrix> bcEqns_reducer;
00400   if (numSlaves > 0) {
00401     bcEqns_reducer.reset(new fei::MatrixReducer(reducer, bcEqns));
00402   }
00403 
00404   fei::Matrix& bcEqns_mat = bcEqns_reducer.get()==NULL ?
00405       *bcEqns : *bcEqns_reducer;
00406 
00407   CHK_ERR( bcManager->finalizeBCEqns(bcEqns_mat, bcs_trump_slaves) );
00408 
00409   if (resolveConflictRequested) {
00410     fei::SharedPtr<fei::FillableMat> mat = bcEqns->getMatrix();
00411     std::vector<int> bcEqnNumbers;
00412     fei::get_row_numbers(*mat, bcEqnNumbers);
00413     CHK_ERR( snl_fei::resolveConflictingCRs(*matrixGraph, bcEqns_mat,
00414                                             bcEqnNumbers) );
00415   }
00416 
00417   std::vector<int> essEqns;
00418   std::vector<double> values;
00419 
00420   std::map<int,fei::FillableMat*>& remotes = bcEqns->getRemotelyOwnedMatrices();
00421   std::map<int,fei::FillableMat*>::iterator
00422     it = remotes.begin(),
00423     it_end = remotes.end();
00424   for(; it!=it_end; ++it) {
00425     fei::impl_utils::separate_BC_eqns( *(it->second), essEqns, values);
00426   }
00427 
00428   CHK_ERR( bcEqns->gatherFromOverlap(false) );
00429 
00430   fei::impl_utils::separate_BC_eqns( *(bcEqns->getMatrix()), essEqns, values);
00431 
00432   if (essEqns.size() > 0) {
00433     int* essEqnsPtr = &essEqns[0];
00434     double* valuesPtr = &values[0];
00435 
00436     for(unsigned i=0; i<essEqns.size(); ++i) {
00437       int eqn = essEqnsPtr[i];
00438       double value = valuesPtr[i];
00439       fei::put_entry(*essBCvalues, eqn, value);
00440     }
00441   }
00442 
00443   return(0);
00444 }
00445 
00446 //----------------------------------------------------------------------------
00447 int snl_fei::LinearSystem_General::implementBCs(bool applyBCs)
00448 {
00449   if (essBCvalues_ != NULL) {
00450     delete essBCvalues_;
00451   }
00452 
00453   essBCvalues_ = new fei::CSVec;
00454 
00455   CHK_ERR( extractDirichletBCs(dbcManager_, matrixGraph_,
00456                        essBCvalues_,  resolveConflictRequested_,
00457                       bcs_trump_slaves_) );
00458 
00459   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00460     FEI_OSTREAM& os = *output_stream_;
00461     std::vector<int>& indices = essBCvalues_->indices();
00462     std::vector<double>& coefs= essBCvalues_->coefs();
00463     for(size_t i=0; i<essBCvalues_->size(); ++i) {
00464       os << "essBCeqns["<<i<<"]: "<<indices[i]<<", "<<coefs[i]<<FEI_ENDL;
00465     }
00466   }
00467 
00468   //If the underlying matrix is a LinearSystemCore instance, then this
00469   //function will return 0, and we're done. A non-zero return-code means
00470   //we should continue and enforce the BCs assuming a general matrix.
00471 
00472   int returncode = enforceEssentialBC_LinSysCore();
00473   if (returncode == 0) {
00474     return(0);
00475   }
00476 
00477   fei::CSVec allEssBCs;
00478   if (!BCenforcement_no_column_mod_) {
00479     fei::impl_utils::global_union(comm_, *essBCvalues_, allEssBCs);
00480 
00481     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00482       FEI_OSTREAM& os = *output_stream_;
00483       os << "  implementBCs, essBCvalues_.size(): "<<essBCvalues_->size()
00484          << ", allEssBCs.size(): " << allEssBCs.size()<<FEI_ENDL;
00485     }
00486   }
00487 
00488   if (essBCvalues_->size() > 0) {
00489     enforceEssentialBC_step_1(*essBCvalues_);
00490   }
00491 
00492   if (!BCenforcement_no_column_mod_ && allEssBCs.size() > 0) {
00493     enforceEssentialBC_step_2(allEssBCs);
00494   }
00495 
00496   return(0);
00497 }
00498 
00499 //----------------------------------------------------------------------------
00500 int snl_fei::LinearSystem_General::enforceEssentialBC_LinSysCore()
00501 {
00502   fei::Matrix* matptr = matrix_.get();
00503   fei::MatrixReducer* matred = dynamic_cast<fei::MatrixReducer*>(matptr);
00504   if (matred != NULL) {
00505     matptr = matred->getTargetMatrix().get();
00506   }
00507 
00508   fei::Matrix_Impl<LinearSystemCore>* lscmatrix =
00509     dynamic_cast<fei::Matrix_Impl<LinearSystemCore>*>(matptr);
00510   if (lscmatrix == 0) {
00511     return(-1);
00512   }
00513 
00514   int localsize = matrixGraph_->getRowSpace()->getNumIndices_Owned();
00515   fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
00516   if (matrixGraph_->getGlobalNumSlaveConstraints() > 0) {
00517     localsize = reducer->getLocalReducedEqns().size();
00518   }
00519 
00520   fei::SharedPtr<fei::FillableMat> inner(new fei::FillableMat);
00521   fei::SharedPtr<fei::Matrix_Impl<fei::FillableMat> > matrix;
00522   matrix.reset(new fei::Matrix_Impl<fei::FillableMat>(inner, matrixGraph_, localsize));
00523 
00524   fei::SharedPtr<fei::SparseRowGraph> remoteGraph =
00525     matrixGraph_->getRemotelyOwnedGraphRows();
00526 
00527   if (!BCenforcement_no_column_mod_) {
00528     CHK_ERR( snl_fei::gatherRemoteEssBCs(*essBCvalues_, remoteGraph.get(), *matrix) );
00529   }
00530 
00531   unsigned numBCRows = inner->getNumRows();
00532 
00533   if (output_stream_ != NULL && output_level_ >= fei::BRIEF_LOGS) {
00534     FEI_OSTREAM& os = *output_stream_;
00535     os << "#enforceEssentialBC_LinSysCore RemEssBCs to enforce: "
00536        << numBCRows << FEI_ENDL;
00537   }
00538 
00539   if (numBCRows > 0 && !BCenforcement_no_column_mod_) {
00540     std::vector<int*> colIndices(numBCRows);
00541     std::vector<double*> coefs(numBCRows);
00542     std::vector<int> colIndLengths(numBCRows);
00543 
00544     fei::CSRMat csrmat(*inner);
00545     fei::SparseRowGraph& srg = csrmat.getGraph();
00546 
00547     int numEqns = csrmat.getNumRows();
00548     int* eqns = &(srg.rowNumbers[0]);
00549     int* rowOffsets = &(srg.rowOffsets[0]);
00550 
00551     for(int i=0; i<numEqns; ++i) {
00552       colIndices[i] = &(srg.packedColumnIndices[rowOffsets[i]]);
00553       coefs[i] = &(csrmat.getPackedCoefs()[rowOffsets[i]]);
00554       colIndLengths[i] = rowOffsets[i+1] - rowOffsets[i];
00555     }
00556 
00557     int** colInds = &colIndices[0];
00558     int* colIndLens = &colIndLengths[0];
00559     double** BCcoefs = &coefs[0];
00560 
00561     if (output_stream_ != NULL && output_level_ > fei::BRIEF_LOGS) {
00562       FEI_OSTREAM& os = *output_stream_;
00563       for(int i=0; i<numEqns; ++i) {
00564         os << "remBCeqn: " << eqns[i] << ", inds/coefs: ";
00565         for(int j=0; j<colIndLens[i]; ++j) {
00566           os << "("<<colInds[i][j]<<","<<BCcoefs[i][j]<<") ";
00567         }
00568         os << FEI_ENDL;
00569       }
00570     }
00571 
00572     int errcode = lscmatrix->getMatrix()->enforceRemoteEssBCs(numEqns,
00573                     eqns,
00574                     colInds,
00575                     colIndLens,
00576                     BCcoefs);
00577     if (errcode != 0) {
00578       return(errcode);
00579     }
00580   }
00581 
00582   int numEqns = essBCvalues_->size();
00583   if (numEqns > 0) {
00584     int* eqns = &(essBCvalues_->indices())[0];
00585     double* bccoefs = &(essBCvalues_->coefs())[0];
00586     std::vector<double> ones(numEqns, 1.0);
00587 
00588     return(lscmatrix->getMatrix()->enforceEssentialBC(eqns, &ones[0],
00589                 bccoefs, numEqns));
00590   }
00591 
00592   return(0);
00593 }
00594 
00595 //----------------------------------------------------------------------------
00596 void snl_fei::LinearSystem_General::enforceEssentialBC_step_1(fei::CSVec& essBCs)
00597 {
00598   //to enforce essential boundary conditions, we do the following:
00599   //
00600   //  1.  for each eqn (== essBCs.indices()[n]), {
00601   //        put zeros in row A[eqn], but leave 1.0 on the diagonal
00602   //        set b[eqn] = essBCs.coefs()[n]
00603   //      }
00604   //
00605   //  2.  for i in 1..numRows (i.e., all rows) {
00606   //        if (i in bcEqns) continue;
00607   //        b[i] -= A[i,eqn] * essBCs.coefs()[n]
00608   //        A[i,eqn] = 0.0;
00609   //      }
00610   //
00611   //It is important to note that for step 1, essBCs need only contain
00612   //local eqns, but for step 2 it should contain *ALL* bc eqns.
00613   //
00614   //This function performs step 1.
00615 
00616   int numEqns = essBCs.size();
00617   int* eqns = &(essBCs.indices())[0];
00618   double* bcCoefs = &(essBCs.coefs())[0];
00619 
00620   std::vector<double> coefs;
00621   std::vector<int> indices;
00622 
00623   fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
00624   bool haveSlaves = reducer.get()!=NULL;
00625 
00626   try {
00627   for(int i=0; i<numEqns; i++) {
00628     int eqn = eqns[i];
00629 
00630     //if slave-constraints are present, the incoming bc-eqns are in
00631     //the reduced equation space. so we actually have to translate them back
00632     //to the unreduced space before passing them into the fei::Matrix object,
00633     //because the fei::Matrix object has machinery to translate unreduced eqns
00634     //to the reduced space.
00635     //Also, our firstLocalOffset_ and lastLocalOffset_ attributes are in the
00636     //unreduced space.
00637     if (haveSlaves) {
00638       eqn = reducer->translateFromReducedEqn(eqn);
00639     }
00640 
00641     if (eqn < firstLocalOffset_ || eqn > lastLocalOffset_) continue;
00642 
00643     //put gamma/alpha on the rhs for this ess-BC equation.
00644     double bcValue = bcCoefs[i];
00645     int err = rhs_->copyIn(1, &eqn, &bcValue);
00646     if (err != 0) {
00647       FEI_OSTRINGSTREAM osstr;
00648       osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_1 ERROR: "
00649       << "err="<<err<<" returned from rhs_->copyIn row="<<eqn;
00650       throw std::runtime_error(osstr.str());
00651     }
00652 
00653     err = getMatrixRow(matrix_.get(), eqn, coefs, indices);
00654     if (err != 0 || indices.size() < 1) {
00655       continue;
00656     }
00657 
00658     int rowLen = indices.size();
00659     int* indPtr = &indices[0];
00660 
00661     //first, put zeros in the row and 1.0 on the diagonal...
00662     for(int j=0; j<rowLen; j++) {
00663       if (indPtr[j] == eqn) coefs[j] = 1.0;
00664       else coefs[j] = 0.0;
00665     }
00666 
00667     double* coefPtr = &coefs[0];
00668 
00669     err = matrix_->copyIn(1, &eqn, rowLen, indPtr, &coefPtr);
00670     if (err != 0) {
00671       FEI_OSTRINGSTREAM osstr;
00672       osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_1 ERROR: "
00673       << "err="<<err<<" returned from matrix_->copyIn row="<<eqn;
00674       throw std::runtime_error(osstr.str());
00675     }
00676   }//for i
00677   }
00678   catch(std::runtime_error& exc) {
00679     FEI_OSTRINGSTREAM osstr;
00680     osstr << "fei::LinearSystem::enforceEssentialBC: ERROR, caught exception: "
00681         << exc.what();
00682     throw std::runtime_error(osstr.str());
00683   }
00684 }
00685 
00686 //----------------------------------------------------------------------------
00687 void snl_fei::LinearSystem_General::enforceEssentialBC_step_2(fei::CSVec& essBCs)
00688 {
00689   //to enforce essential boundary conditions, we do the following:
00690   //
00691   //  1.  for each eqn (== essBCs.indices()[n]), {
00692   //        put zeros in row A[eqn], but leave 1.0 on the diagonal
00693   //        set b[eqn] = essBCs.coefs()[n]
00694   //      }
00695   //
00696   //  2.  for i in 1..numRows (i.e., all rows) {
00697   //        if (i in bcEqns) continue;
00698   //        b[i] -= A[i,eqn] * essBCs.coefs()[n]
00699   //        A[i,eqn] = 0.0;
00700   //      }
00701   //
00702   //It is important to note that for step 1, essBCs need only contain
00703   //local eqns, but for step 2 it should contain *ALL* bc eqns.
00704   //
00705   //This function performs step 2.
00706 
00707   int numBCeqns = essBCs.size();
00708   if (numBCeqns < 1) {
00709     return;
00710   }
00711 
00712   int* bcEqns = &(essBCs.indices())[0];
00713   double* bcCoefs = &(essBCs.coefs())[0];
00714 
00715   fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
00716   bool haveSlaves = reducer.get()!=NULL;
00717   if (haveSlaves) {
00718     for(int i=0; i<numBCeqns; ++i) {
00719       bcEqns[i] = reducer->translateFromReducedEqn(bcEqns[i]);
00720     }
00721   }
00722 
00723   int firstBCeqn = bcEqns[0];
00724   int lastBCeqn = bcEqns[numBCeqns-1];
00725 
00726   std::vector<double> coefs;
00727   std::vector<int> indices;
00728 
00729   int insertPoint;
00730 
00731   int nextBCeqnOffset = 0;
00732   int nextBCeqn = bcEqns[nextBCeqnOffset];
00733 
00734   for(int i=firstLocalOffset_; i<=lastLocalOffset_; ++i) {
00735     if (haveSlaves) {
00736       if (reducer->isSlaveEqn(i)) continue;
00737     }
00738 
00739     bool should_continue = false;
00740     if (i >= nextBCeqn) {
00741       if (i == nextBCeqn) {
00742   ++nextBCeqnOffset;
00743   if (nextBCeqnOffset < numBCeqns) {
00744     nextBCeqn = bcEqns[nextBCeqnOffset];
00745   }
00746   else {
00747     nextBCeqn = lastLocalOffset_+1;
00748   }
00749 
00750   should_continue = true;
00751       }
00752       else {
00753   while(nextBCeqn <= i) {
00754     if (nextBCeqn == i) should_continue = true;
00755     ++nextBCeqnOffset;
00756     if (nextBCeqnOffset < numBCeqns) {
00757       nextBCeqn = bcEqns[nextBCeqnOffset];
00758     }
00759     else {
00760       nextBCeqn = lastLocalOffset_+1;
00761     }
00762   }
00763       }
00764     }
00765 
00766     if (should_continue) continue;
00767 
00768     int err = getMatrixRow(matrix_.get(), i, coefs, indices);
00769     if (err != 0 || indices.size() < 1) {
00770       continue;
00771     }
00772 
00773     int numIndices = indices.size();
00774     int* indicesPtr = &indices[0];
00775     double* coefsPtr = &coefs[0];
00776     bool modifiedCoef = false;
00777 
00778     fei::insertion_sort_with_companions(numIndices, indicesPtr, coefsPtr);
00779 
00780     if (indicesPtr[0] > lastBCeqn || indicesPtr[numIndices-1] < firstBCeqn) {
00781       continue;
00782     }
00783 
00784     double value = 0.0;
00785     int offset = 0;
00786 
00787     for(int j=0; j<numIndices; ++j) {
00788       int idx = indicesPtr[j];
00789       offset = fei::binarySearch(idx, bcEqns, numBCeqns,
00790              insertPoint);
00791       if (offset > -1) {
00792   value -= bcCoefs[offset]*coefsPtr[j];
00793 
00794   coefsPtr[j] = 0.0;
00795   modifiedCoef = true;
00796       }
00797     }
00798 
00799     if (modifiedCoef) {
00800       err = matrix_->copyIn(1, &i, numIndices, indicesPtr, &coefsPtr);
00801       if (err != 0) {
00802   FEI_OSTRINGSTREAM osstr;
00803   osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_2 ERROR: "
00804         << "err="<<err<<" returned from matrix_->copyIn, row="<<i;
00805   throw std::runtime_error(osstr.str());
00806       }
00807     }
00808 
00809     const double fei_eps = 1.e-49;
00810     if (std::abs(value) > fei_eps) {
00811       rhs_->sumIn(1, &i, &value);
00812 
00813       if (output_level_ >= fei::FULL_LOGS && output_stream_ != 0) {
00814   FEI_OSTREAM& os = *output_stream_;
00815   os << "enfEssBC_step2: rhs["<<i<<"] += "<<value<<FEI_ENDL;
00816       }
00817     }
00818   }
00819 }
00820 
00821 //----------------------------------------------------------------------------
00822 int snl_fei::LinearSystem_General::getMatrixRow(fei::Matrix* matrix, int row,
00823             std::vector<double>& coefs,
00824             std::vector<int>& indices)
00825 {
00826   int len = 0;
00827   int err = matrix->getRowLength(row, len);
00828   if (err != 0) {
00829     coefs.resize(0);
00830     indices.resize(0);
00831     return(err);
00832   }
00833 
00834   if ((int)coefs.size() != len) {
00835     coefs.resize(len);
00836   }
00837   if ((int)indices.size() != len) {
00838     indices.resize(len);
00839   }
00840 
00841   CHK_ERR( matrix->copyOutRow(row, len, &coefs[0], &indices[0]));
00842 
00843   return(0);
00844 }
00845 
00846 //----------------------------------------------------------------------------
00847 int snl_fei::LinearSystem_General::loadLagrangeConstraint(int constraintID,
00848                 const double *weights,
00849                 double rhsValue)
00850 {
00851   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00852     FEI_OSTREAM& os = *output_stream_;
00853     os << "loadLagrangeConstraint crID: "<<constraintID<<FEI_ENDL;
00854   }
00855 
00856   Constraint<fei::Record<int>*>* cr =
00857     matrixGraph_->getLagrangeConstraint(constraintID);
00858   if (cr == NULL) {
00859     return(-1);
00860   }
00861 
00862   CHK_ERR( matrixGraph_->getConstraintConnectivityIndices(cr, iwork_) );
00863 
00864   //Let's attach the weights to the constraint-record now.
00865   std::vector<double>& cr_weights = cr->getMasterWeights();
00866   cr_weights.resize(iwork_.size());
00867   for(unsigned i=0; i<iwork_.size(); ++i) {
00868     cr_weights.push_back(weights[i]);
00869   }
00870 
00871   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph_->getRowSpace();
00872 
00873   int crEqn = -1;
00874   CHK_ERR( vecSpace->getGlobalIndex(cr->getIDType(),
00875              cr->getConstraintID(),
00876              crEqn) );
00877 
00878   //now add the row contribution to the matrix and rhs
00879   int numIndices = iwork_.size();
00880   int* indicesPtr = &(iwork_[0]);
00881 
00882   CHK_ERR( matrix_->sumIn(1, &crEqn, numIndices, indicesPtr, &weights) );
00883 
00884   CHK_ERR( rhs_->sumIn(1, &crEqn, &rhsValue) );
00885 
00886   //now add the column contributions to the matrix
00887   for(int k=0; k<numIndices; ++k) {
00888     double* thisWeight = (double*)(&(weights[k]));
00889     CHK_ERR( matrix_->sumIn(1, &(indicesPtr[k]), 1, &crEqn, &thisWeight) );
00890   }
00891 
00892   return(0);
00893 }
00894 
00895 //----------------------------------------------------------------------------
00896 int snl_fei::LinearSystem_General::loadPenaltyConstraint(int constraintID,
00897                const double *weights,
00898                double penaltyValue,
00899                double rhsValue)
00900 {
00901   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00902     FEI_OSTREAM& os = *output_stream_;
00903     os << "loadPenaltyConstraint crID: "<<constraintID<<FEI_ENDL;
00904   }
00905 
00906   Constraint<fei::Record<int>*>* cr =
00907     matrixGraph_->getPenaltyConstraint(constraintID);
00908   if (cr == NULL) {
00909     return(-1);
00910   }
00911 
00912   CHK_ERR( matrixGraph_->getConstraintConnectivityIndices(cr, iwork_) );
00913 
00914   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph_->getRowSpace();
00915 
00916   int numIndices = iwork_.size();
00917   int* indicesPtr = &(iwork_[0]);
00918 
00919   //now add the contributions to the matrix and rhs
00920   std::vector<double> coefs(numIndices);
00921   double* coefPtr = &coefs[0];
00922   for(int i=0; i<numIndices; ++i) {
00923     for(int j=0; j<numIndices; ++j) {
00924       coefPtr[j] = weights[i]*weights[j]*penaltyValue;
00925     }
00926     CHK_ERR( matrix_->sumIn(1, &(indicesPtr[i]), numIndices, indicesPtr,
00927           &coefPtr) );
00928 
00929     double rhsCoef = weights[i]*penaltyValue*rhsValue;
00930     CHK_ERR( rhs_->sumIn(1, &(indicesPtr[i]), &rhsCoef) );
00931   }
00932 
00933   return(0);
00934 }
00935 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends