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