snl_fei_Utils.cpp

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

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