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