FEI Version of the Day
snl_fei_Utils.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 
00046 #include <limits>
00047 #include <cmath>
00048 #include <stdexcept>
00049 
00050 #include <snl_fei_Utils.hpp>
00051 #include "fei_Record.hpp"
00052 #include <fei_MatrixGraph.hpp>
00053 #include <fei_SparseRowGraph.hpp>
00054 #include <fei_Matrix_Impl.hpp>
00055 #include <fei_ParameterSet.hpp>
00056 
00057 #include <fei_CommUtils.hpp>
00058 #include <fei_TemplateUtils.hpp>
00059 #include <fei_CSVec.hpp>
00060 #include <fei_chk_mpi.hpp>
00061 
00062 #undef fei_file
00063 #define fei_file "snl_fei_Utils.cpp"
00064 #include <fei_ErrMacros.hpp>
00065 
00066 //----------------------------------------------------------------------------
00067 const char* snl_fei::getParam(const char* key,
00068             int numParams,
00069             const char* const* paramStrings)
00070 {
00071   int foundOffset = -1;
00072   return( getParam(key, numParams, paramStrings, foundOffset) );
00073 }
00074 
00075 //----------------------------------------------------------------------------
00076 const char* snl_fei::getParamValue(const char* key,
00077            int numParams,
00078            const char* const* paramStrings,
00079            char separator)
00080 {
00081   const char* param = getParam(key, numParams, paramStrings);
00082 
00083   return( param != NULL ? skipSeparator(param, separator) : NULL );
00084 }
00085 
00086 //----------------------------------------------------------------------------
00087 const char* snl_fei::getParamValue(const char* key,
00088            int numParams,
00089            const char* const* paramStrings,
00090            int& foundOffset,
00091            char separator)
00092 {
00093   const char* param = getParam(key, numParams, paramStrings, foundOffset);
00094 
00095   return( param != NULL ? skipSeparator(param, separator) : NULL );
00096 }
00097 
00098 //----------------------------------------------------------------------------
00099 int snl_fei::getIntParamValue(const char* key,
00100             int numParams,
00101             const char*const* params,
00102             int& paramValue)
00103 {
00104   const char* tempstr = getParamValue(key, numParams, params);
00105   if (tempstr==NULL) return(-1);
00106 
00107   std::string strg(tempstr);
00108   FEI_ISTRINGSTREAM isstr(strg);
00109   isstr >> paramValue;
00110   return( isstr.fail() ? -1 : 0 );
00111 }
00112 
00113 //----------------------------------------------------------------------------
00114 int snl_fei::getDoubleParamValue(const char* key,
00115          int numParams,
00116          const char*const* params,
00117          double& paramValue)
00118 {
00119   const char* tempstr = getParamValue(key, numParams, params);
00120   if (tempstr==NULL) return(-1);
00121 
00122   std::string strg(tempstr);
00123   FEI_ISTRINGSTREAM isstr(strg);
00124   isstr >> paramValue;
00125   return( isstr.fail() ? -1 : 0 );
00126 }
00127 
00128 //----------------------------------------------------------------------------
00129 int snl_fei::getDoubleParamValue(const char* key,
00130          std::vector<std::string>& params,
00131          double& paramValue)
00132 {
00133   const char* tempstr = getParamValue(key, params);
00134   if (tempstr==NULL) return(-1);
00135 
00136   std::string strg(tempstr);
00137   FEI_ISTRINGSTREAM isstr(strg);
00138   isstr >> paramValue;
00139   return( isstr.fail() ? -1 : 0 );
00140 }
00141 
00142 //----------------------------------------------------------------------------
00143 const char* snl_fei::getParam(const char* key,
00144             int numParams,
00145             const char* const* paramStrings,
00146             int& foundOffset)
00147 {
00148   const char* returnPtr = NULL;
00149   foundOffset = -1;
00150 
00151   //if the search-key is null, or if the list of strings to be searched is null,
00152   //or if the number of strings to be searched is null, then return now.
00153   //
00154   if (key == NULL || paramStrings == NULL || numParams == 0) {
00155     return(returnPtr);
00156   }
00157 
00158   unsigned keyLen = strlen(key);
00159 
00160   for(int i=numParams-1; i>=0; --i) {
00161     const char* paramStr = paramStrings[i];
00162 
00163     //if the i-th paramString is null, skip it.
00164     if (paramStr == NULL) continue;
00165 
00166     unsigned paramStr_len = leading_substring_length(paramStr);
00167     if (paramStr_len != keyLen) continue;
00168 
00169     if (strncmp(key, paramStr, keyLen) == 0) {
00170       returnPtr = paramStr;
00171       foundOffset = i;
00172       return(returnPtr);
00173     }
00174   }
00175 
00176   return(returnPtr);
00177 }
00178 
00179 //----------------------------------------------------------------------------
00180 const char* snl_fei::getParam(const char* key,
00181             std::vector<std::string>& params,
00182             int& foundOffset)
00183 {
00184   std::vector<std::string>::iterator
00185     p_iter = params.begin(),
00186     p_end = params.end();
00187 
00188   int offset = 0;
00189   for(; p_iter != p_end; ++p_iter, ++offset) {
00190     std::string& pstr = *p_iter;
00191     int ssize = pstr.size();
00192     if (ssize < 1) continue;
00193 
00194     std::string::size_type i = pstr.find(key);
00195     if (i == 0) {
00196       foundOffset = offset;
00197       return( pstr.c_str() );
00198     }
00199   }
00200 
00201   return( NULL );
00202 }
00203 
00204 //----------------------------------------------------------------------------
00205 const char* snl_fei::getParamValue(const char* key,
00206            std::vector<std::string>& params,
00207            char separator)
00208 {
00209   int offset = 0;
00210   return( skipSeparator( getParam(key, params, offset), separator) );
00211 }
00212 
00213 //----------------------------------------------------------------------------
00214 void snl_fei::separate_string(const char* input_string,
00215             const char* substring,
00216             const char*& before_substring,
00217             int& len_before_substring,
00218             const char*& after_substring,
00219             int& len_after_substring)
00220 {
00221   if (input_string == NULL) {
00222     before_substring = NULL;
00223     len_before_substring = 0;
00224     after_substring = NULL;
00225     len_after_substring = 0;
00226     return;
00227   }
00228 
00229   std::string std_input_string(input_string);
00230   int len_input_string = std_input_string.size();
00231   before_substring = input_string;
00232 
00233   if (substring == NULL) {
00234     len_before_substring = len_input_string;
00235     after_substring = NULL;
00236     len_after_substring = 0;
00237     return;
00238   }
00239 
00240   std::string std_substring(substring);
00241   int len_substring = std_substring.size();
00242 
00243   size_t found = std_input_string.find(std_substring);
00244   if (found == std::string::npos) {
00245     len_before_substring = len_input_string;
00246     after_substring = NULL;
00247     len_after_substring = 0;
00248     return;
00249   }
00250 
00251   size_t rfound = std_input_string.rfind(std_substring);
00252   after_substring = input_string+rfound+std_substring.size();
00253   len_before_substring = found;
00254   len_after_substring = len_input_string - rfound - len_substring;
00255 }
00256 
00257 //----------------------------------------------------------------------------
00258 int snl_fei::storeNamedAttribute(const char* name,
00259          void* attribute,
00260          std::vector<char*>& attributeNames,
00261          std::vector<void*>& attributes)
00262 {
00263   int offset = -1;
00264   const char* foundname = attributeNames.empty() ?
00265     0 : snl_fei::getParam(name, attributeNames.size(),
00266         &(attributeNames[0]), offset);
00267 
00268   if (foundname != 0) {
00269     attributes[offset] = attribute;
00270   }
00271   else {
00272     char* newname = new char[strlen(name)+1];
00273     strcpy(newname, name);
00274     attributeNames.push_back(newname);
00275     attributes.push_back(attribute);
00276   }
00277 
00278   return(0);
00279 }
00280 
00281 //----------------------------------------------------------------------------
00282 void* snl_fei::retrieveNamedAttribute(const char* name,
00283               std::vector<char*>& attributeNames,
00284               std::vector<void*>& attributes)
00285 {
00286   int offset = -1;
00287   const char* foundname = attributeNames.empty() ?
00288     0 : snl_fei::getParam(name, attributeNames.size(),
00289         &(attributeNames[0]), offset);
00290 
00291   void* returnVal = NULL;
00292 
00293   if (foundname != 0) {
00294     returnVal = attributes[offset];
00295   }
00296   else {
00297     return( returnVal );
00298   }
00299 
00300   return( returnVal );
00301 }
00302 
00303 //----------------------------------------------------------------------------
00304 unsigned snl_fei::leading_substring_length(const char* string)
00305 {
00306   if (string == NULL) {
00307     return(0);
00308   }
00309 
00310   const char* lastchar = string;
00311   unsigned len = 0;
00312 
00313   while(*lastchar != ' ' && *lastchar != '\t' && *lastchar != '\0') {
00314     ++len;
00315     ++lastchar;
00316   }
00317 
00318   return(len);
00319 }
00320 
00321 //----------------------------------------------------------------------------
00322 const char* snl_fei::skipSeparator(const char* paramString,
00323            char separator)
00324 {
00325   if (paramString == NULL) return(NULL);
00326 
00327   const char* result = strchr(paramString, separator);
00328 
00329   if (result == NULL) return(result);
00330 
00331   //allow for the possibility that separator is repeated
00332   while(result[0] == separator) result++;
00333 
00334   return( result );
00335 }
00336 
00337 //----------------------------------------------------------------------------
00338 int snl_fei::mergeStringLists(char**& strings,
00339             int& numStrings,
00340             const char*const* stringsToMerge,
00341             int numStringsToMerge)
00342 {
00343   int i;
00344   if (numStrings == 0) {
00345     strings = new char*[numStringsToMerge];
00346 
00347     for(i=0; i<numStringsToMerge; ++i){
00348       strings[i] = new char[strlen(stringsToMerge[i])+1];
00349 
00350       strcpy(strings[i], stringsToMerge[i]);
00351     }
00352 
00353     numStrings = numStringsToMerge;
00354   }
00355   else {
00356     int numMerged = 0;
00357     bool* merged = new bool[numStringsToMerge];
00358     for(i=0; i<numStringsToMerge; ++i) {
00359       const char* temp;
00360       const char* temp2;
00361       int len_temp, len_temp2;
00362       separate_string(stringsToMerge[i], " ",
00363                       temp, len_temp, temp2, len_temp2);
00364       std::string strtemp(temp, len_temp);
00365       int foundOffset = -1;
00366       const char* matchingString = 
00367   snl_fei::getParam(strtemp.c_str(), numStrings, strings, foundOffset);
00368       if (matchingString != NULL) {
00369         int len = strlen(stringsToMerge[i])+1;
00370   delete [] strings[foundOffset];
00371   strings[foundOffset] = new char[len];
00372   strcpy(strings[foundOffset], stringsToMerge[i]);
00373   merged[i] = true;
00374   ++numMerged;
00375       }
00376       else merged[i] = false;
00377     }
00378 
00379     if (numMerged == numStringsToMerge) {
00380       delete [] merged;
00381       return(0);
00382     }
00383 
00384     char** newStrings = new char*[numStrings+numStringsToMerge-numMerged];
00385     for(i=0; i<numStrings; ++i) {
00386       newStrings[i] = strings[i];
00387     }
00388     int offset = numStrings;
00389     for(i=0; i<numStringsToMerge; ++i) {
00390       if (!merged[i]) {
00391   int len = strlen(stringsToMerge[i])+1;
00392   newStrings[offset] = new char[len];
00393   strcpy(newStrings[offset++], stringsToMerge[i]);
00394       }
00395     }
00396     delete [] merged;
00397 
00398     delete [] strings;
00399     strings = newStrings;
00400     numStrings+= numStringsToMerge-numMerged;
00401   }
00402 
00403   return(0);
00404 }
00405 
00406 //----------------------------------------------------------------------------
00407 int snl_fei::resolveConflictingCRs(fei::MatrixGraph& matrixGraph,
00408            fei::Matrix& bcEqns,
00409                                    const std::vector<int>& bcEqnNumbers)
00410 {
00411   int numLagrangeConstraints = matrixGraph.getLocalNumLagrangeConstraints();
00412   if (numLagrangeConstraints < 1) {
00413     return(0);
00414   }
00415 
00416   double coefs[3];
00417   int indices[3];
00418   indices[0] = 0;
00419   indices[1] = 1;
00420   indices[2] = 2;
00421   coefs[0] = 1.0;
00422   coefs[1] = 0.0;
00423   coefs[2] = 1.0;
00424 
00425   double* coefPtr = &(coefs[0]);
00426   int* indicesPtr = &(indices[0]);
00427 
00428   double fei_eps = 1.e-49;
00429 
00430   std::vector<int> cr_indices;
00431   std::map<int,Constraint<fei::Record<int>*>*>& lagrangeConstraints =
00432     matrixGraph.getLagrangeConstraints();
00433 
00434   std::map<int,Constraint<fei::Record<int>*>*>::const_iterator
00435     cr_iter = lagrangeConstraints.begin(),
00436     cr_end  = lagrangeConstraints.end();
00437 
00438   while(cr_iter != cr_end) {
00439     Constraint<fei::Record<int>*>* cr = (*cr_iter).second;
00440 
00441     CHK_ERR( matrixGraph.getConstraintConnectivityIndices(cr, cr_indices) );
00442 
00443     std::vector<double>& weights = cr->getMasterWeights();
00444     double* weightsPtr = &weights[0];
00445 
00446     int len = cr_indices.size();
00447     int* cr_indPtr = &(cr_indices[0]);
00448     for(int j=0; j<len; ++j) {
00449       if (std::abs(weightsPtr[j] + 1.0) > fei_eps) continue;
00450 
00451       int eqn = cr_indPtr[j];
00452       if (fei::binarySearch(eqn, bcEqnNumbers) > -1) {
00453   int cr_eqn = cr->getEqnNumber();
00454 
00455   CHK_ERR( bcEqns.copyIn(1, &cr_eqn, 3, indicesPtr,
00456              (double**)(&coefPtr)) );
00457       }
00458     }
00459     ++cr_iter;
00460   }
00461 
00462   return(0);
00463 }
00464 
00465 //----------------------------------------------------------------------------
00466 int snl_fei::gatherRemoteEssBCs(fei::CSVec& essBCs,
00467         fei::SparseRowGraph* remoteGraph,
00468         fei::Matrix& matrix)
00469 {
00470   std::vector<int>& rrowOffs = remoteGraph->rowOffsets;
00471   std::vector<int>& rcols = remoteGraph->packedColumnIndices;
00472 
00473   int numEssEqns = essBCs.size();
00474   if (numEssEqns > 0) {
00475     int* essEqns = &(essBCs.indices()[0]);
00476     double* coefs = &(essBCs.coefs()[0]);
00477 
00478     if (rrowOffs.size() > 0 && rcols.size() > 0) {
00479 
00480       int* rowOffsPtr = &(rrowOffs[0]);
00481       int* rcolsPtr = &(rcols[0]);
00482 
00483       for(int j=0; j<numEssEqns; ++j) {
00484 
00485         int eqn = essEqns[j];
00486 
00487         for(unsigned i=0; i<rrowOffs.size()-1; ++i) {
00488           int len = rowOffsPtr[i+1]-rowOffsPtr[i];
00489           int* colsPtr = &(rcolsPtr[rowOffsPtr[i]]);
00490 
00491           if (fei::binarySearch(eqn, colsPtr, len) > -1) {
00492             double coef = coefs[j];
00493             double* coefPtr = &coef;
00494 
00495             for(int k=0; k<len; ++k) {
00496               CHK_ERR( matrix.copyIn(1, &(colsPtr[k]),
00497                        1, &(eqn), &coefPtr) );
00498             }
00499           }
00500         }
00501       }
00502     }
00503   }
00504 
00505   CHK_ERR( matrix.gatherFromOverlap(false) );
00506 
00507   return(0);
00508 }
00509 
00510 fei::SharedPtr<fei::SparseRowGraph>
00511 snl_fei::mergeSparseRowGraphs(const fei::SparseRowGraph* srg1,
00512                               const fei::SparseRowGraph* srg2)
00513 {
00514   fei::SharedPtr<fei::SparseRowGraph> newgraph(new fei::SparseRowGraph);
00515 
00516   int numrows = srg1->rowNumbers.size() + srg2->rowNumbers.size();
00517   int nnz = srg1->packedColumnIndices.size() + srg2->packedColumnIndices.size();
00518 
00519   newgraph->rowNumbers.resize(numrows);
00520   newgraph->rowOffsets.resize(numrows+1);
00521   newgraph->packedColumnIndices.resize(nnz);
00522 
00523   std::map<int,int> rowmap;
00524 
00525   for(unsigned i=0; i<srg1->rowNumbers.size(); ++i) {
00526     int rowlen = srg1->rowOffsets[i+1]-srg1->rowOffsets[i];
00527     rowmap.insert(std::make_pair(srg1->rowNumbers[i], rowlen));
00528   }
00529 
00530   for(unsigned i=0; i<srg2->rowNumbers.size(); ++i) {
00531     int rowlen = srg2->rowOffsets[i+1]-srg2->rowOffsets[i];
00532     rowmap.insert(std::make_pair(srg2->rowNumbers[i], rowlen));
00533   }
00534 
00535   std::map<int,int>::iterator
00536     r_iter = rowmap.begin(),
00537     r_end  = rowmap.end();
00538 
00539   int offset = 0;
00540   for(unsigned i=0; r_iter != r_end; ++r_iter, ++i) {
00541     newgraph->rowNumbers[i] = r_iter->first;
00542     int rowlen = r_iter->second;
00543     newgraph->rowOffsets[i] = offset;
00544     offset += rowlen;
00545     r_iter->second = i;
00546   }
00547   newgraph->rowOffsets[numrows] = offset;
00548 
00549   for(unsigned i=0; i<srg1->rowNumbers.size(); ++i) {
00550     r_iter = rowmap.find(srg1->rowNumbers[i]);
00551 
00552     int newcoloffset = newgraph->rowOffsets[r_iter->second];
00553 
00554     int jbegin = srg1->rowOffsets[i];
00555     int jend   = srg1->rowOffsets[i+1];
00556     for(int j=jbegin; j<jend; ++j) {
00557       newgraph->packedColumnIndices[newcoloffset++] =
00558         srg1->packedColumnIndices[j];
00559     }
00560   }
00561 
00562   for(unsigned i=0; i<srg2->rowNumbers.size(); ++i) {
00563     r_iter = rowmap.find(srg2->rowNumbers[i]);
00564 
00565     int newcoloffset = newgraph->rowOffsets[r_iter->second];
00566 
00567     int jbegin = srg2->rowOffsets[i];
00568     int jend   = srg2->rowOffsets[i+1];
00569     for(int j=jbegin; j<jend; ++j) {
00570       newgraph->packedColumnIndices[newcoloffset++] =
00571         srg2->packedColumnIndices[j];
00572     }
00573   }
00574 
00575   return(newgraph);
00576 }
00577 
00578 void snl_fei::copy2DBlockDiagToColumnContig(int numBlocks,
00579               const int* blockSizes,
00580               const double*const* values2d,
00581               int format,
00582               double* colcontigvalues)
00583 {
00584   int i, j, k, offset, coffset;
00585 
00586   switch(format) {
00587   case FEI_BLOCK_DIAGONAL_ROW:
00588     offset = 0;
00589     coffset = 0;
00590     for(i=0; i<numBlocks; ++i) {
00591       int bsize = blockSizes[i];
00592       for(j=0; j<bsize; ++j) {
00593   for(k=0; k<bsize; ++k) {
00594     colcontigvalues[coffset+j*bsize+k] = values2d[offset+j][k];
00595   }
00596       }
00597       offset += bsize;
00598       coffset += bsize*bsize;
00599     }
00600     break;
00601   default:
00602     FEI_OSTRINGSTREAM osstr;
00603     osstr << "copy2DBlockDiagToColumnContig ERROR, format="<<format
00604     << " not recognized.";
00605     throw std::runtime_error(osstr.str());
00606   }
00607 }
00608 
00609 void snl_fei::copy2DToColumnContig(int numrows,
00610            int numcols,
00611            const double*const* values2d,
00612            int format,
00613            double* colcontigvalues)
00614 {
00615   int i, j;
00616 
00617   switch(format) {
00618 
00619   case FEI_DENSE_ROW:
00620     for(i=0; i<numrows; ++i) {
00621       for(j=0; j<numcols; ++j) {
00622   colcontigvalues[j*numrows+i] = values2d[i][j];
00623       }
00624     }
00625     break;
00626 
00627   case FEI_DENSE_COL:
00628     for(j=0; j<numcols; ++j) {
00629       for(i=0; i<numrows; ++i) {
00630   colcontigvalues[j*numrows+i] = values2d[j][i];
00631       }
00632     }
00633     break;
00634 
00635   default:
00636     FEI_OSTRINGSTREAM osstr;
00637     osstr << "copy2DToColumnContig ERROR, format="<<format<<" not recognized.";
00638     throw std::runtime_error(osstr.str());
00639   }
00640 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends