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

Generated on Tue Jul 13 09:27:46 2010 for FEI by  doxygen 1.4.7