fei_impl_utils.cpp

00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2008 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_CommUtils.hpp>
00010 #include <fei_iostream.hpp>
00011 #include <fei_impl_utils.hpp>
00012 #include <fei_FillableVec.hpp>
00013 #include <fei_FillableMat.hpp>
00014 #include <fei_CSRMat.hpp>
00015 #include <fei_CSVec.hpp>
00016 #include <fei_Graph.hpp>
00017 #include <fei_Matrix.hpp>
00018 #include <fei_Reducer.hpp>
00019 
00020 namespace fei {
00021 namespace impl_utils {
00022 
00023 //----------------------------------------------------------------------------
00024 void find_offsets(const std::vector<int>& sources,
00025                   const std::vector<int>& targets,
00026                   std::vector<int>& offsets)
00027 {
00028   offsets.assign(sources.size(), -1);
00029 
00030   size_t ssize = sources.size(), tsize = targets.size();
00031   size_t si = 0, ti = 0;
00032 
00033   const int* sourcesptr = &sources[0];
00034   const int* targetsptr = &targets[0];
00035 
00036   while(si<ssize && ti<tsize) {
00037     if (sourcesptr[si] == targetsptr[ti]) {
00038       offsets[si++] = ti++;
00039       continue;
00040     }
00041 
00042     while(sourcesptr[si] < targetsptr[ti] && si<ssize) {
00043       ++si;
00044     }
00045 
00046     while(sourcesptr[si] > targetsptr[ti] && ti<tsize) {
00047       ++ti;
00048     }
00049   }
00050 }
00051 
00052 //----------------------------------------------------------------------------
00053 void pack_FillableMat(const fei::FillableMat& mat,
00054                       std::vector<int>& intdata,
00055                       std::vector<double>& doubledata)
00056 {
00057   int nrows = mat.getNumRows();
00058   int nnz = fei::count_nnz(mat);
00059 
00060   intdata.resize(1 + nrows*2 + nnz);
00061   doubledata.resize(nnz);
00062 
00063   int ioffset = 0;
00064   int doffset = 0;
00065 
00066   intdata[ioffset++] = nrows;
00067 
00068   fei::FillableMat::const_iterator
00069     r_iter = mat.begin(),
00070     r_end = mat.end();
00071 
00072   for(; r_iter!=r_end; ++r_iter) {
00073     int rowNumber = r_iter->first;
00074     const fei::FillableVec* row = r_iter->second;
00075 
00076     intdata[ioffset++] = rowNumber;
00077     intdata[ioffset++] = row->size();
00078 
00079     fei::FillableVec::const_iterator
00080       iter = row->begin(),
00081       iend = row->end();
00082     for(; iter!=iend; ++iter) {
00083       intdata[ioffset++] = iter->first;
00084       doubledata[doffset++] = iter->second;
00085     }
00086   }
00087 }
00088 
00089 //----------------------------------------------------------------------------
00090 void unpack_FillableMat(const std::vector<int>& intdata,
00091                         const std::vector<double>& doubledata,
00092                         fei::FillableMat& mat,
00093                         bool clear_mat_on_entry,
00094                         bool overwrite_entries)
00095 {
00096   if (clear_mat_on_entry) {
00097     mat.clear();
00098   }
00099 
00100   if (intdata.size() < 1) {
00101     return;
00102   }
00103 
00104   int ioffset = 0;
00105   int doffset = 0;
00106 
00107   int nrows = intdata[ioffset++];
00108 
00109   for(int i=0; i<nrows; ++i) {
00110     int row = intdata[ioffset++];
00111     int rowlen = intdata[ioffset++];
00112 
00113     for(int j=0; j<rowlen; ++j) {
00114       int col = intdata[ioffset++];
00115       double coef = doubledata[doffset++];
00116 
00117       if (overwrite_entries) {
00118         mat.putCoef(row, col, coef);
00119       }
00120       else {
00121         mat.sumInCoef(row, col, coef);
00122       }
00123     }
00124   }
00125 }
00126 
00127 //----------------------------------------------------------------------------
00128 void separate_BC_eqns(const fei::FillableMat& mat,
00129                     std::vector<int>& bcEqns,
00130                     std::vector<double>& bcVals)
00131 {
00132   fei::FillableMat::const_iterator
00133     m_iter = mat.begin(),
00134     m_end = mat.end();
00135 
00136   for(; m_iter != m_end; ++m_iter) {
00137     int eqn = m_iter->first;
00138     const fei::FillableVec* row = m_iter->second;
00139 
00140     std::vector<int>::iterator
00141       iter = std::lower_bound(bcEqns.begin(), bcEqns.end(), eqn);
00142     if (iter == bcEqns.end() || *iter != eqn) {
00143       size_t offset = iter - bcEqns.begin();
00144       bcEqns.insert(iter, eqn);
00145       std::vector<double>::iterator viter = bcVals.begin();
00146       viter += offset;
00147       try {
00148         double val = row->getEntry(eqn);
00149         bcVals.insert(viter, val);
00150       }
00151       catch(...) {
00152         FEI_CERR << "separate_BC_eqns ERROR, failed to find coef for eqn " << eqn;
00153         throw;
00154       }
00155     }
00156   }
00157 }
00158 
00159 //----------------------------------------------------------------------------
00160 void create_col_to_row_map(const fei::FillableMat& mat,
00161                            std::multimap<int,int>& crmap)
00162 {
00163   crmap.clear();
00164 
00165   if (mat.getNumRows() == 0) return;
00166 
00167   std::vector<int> rowNumbers;
00168   fei::get_row_numbers(mat, rowNumbers);
00169 
00170   fei::FillableMat::const_iterator
00171     m_iter = mat.begin(),
00172     m_end = mat.end();
00173 
00174   for(; m_iter != m_end; ++m_iter) {
00175     int rowNum = m_iter->first;
00176     const fei::FillableVec* rowvec = m_iter->second;
00177 
00178     fei::FillableVec::const_iterator
00179       r_iter = rowvec->begin(),
00180       r_end = rowvec->end();
00181 
00182     for(; r_iter != r_end; ++r_iter) {
00183       int colNum = r_iter->first;
00184 
00185       crmap.insert(std::make_pair(colNum, rowNum));
00186     }
00187   }
00188 }
00189 
00190 //----------------------------------------------------------------------------
00191 int remove_couplings(fei::FillableMat& mat)
00192 {
00193   int levelsOfCoupling = 0;
00194 
00195   std::vector<int> rowNumbers;
00196   fei::get_row_numbers(mat, rowNumbers);
00197 
00198   bool finished = false;
00199   while(!finished) {
00200     std::multimap<int,int> crmap;
00201     create_col_to_row_map(mat, crmap);
00202 
00203     typedef std::multimap<int,int>::iterator MM_Iter;
00204 
00205     fei::FillableMat::iterator
00206       m_iter = mat.begin(),
00207       m_end = mat.end();
00208 
00209     bool foundCoupling = false;
00210     for(; m_iter != m_end; ++m_iter) {
00211       int rownum = m_iter->first;
00212       fei::FillableVec* mrow = m_iter->second;
00213 
00214       //now find which rows contain 'rownum' as a column-index:
00215       std::pair<MM_Iter,MM_Iter> mmi = crmap.equal_range(rownum);
00216 
00217       //loop over the rows which contain 'rownum' as a column-index:
00218       for(MM_Iter cri = mmi.first; cri != mmi.second; ++cri) {
00219         int cri_row = cri->second;
00220 
00221         fei::FillableVec* frow = mat.getRow(cri_row);
00222 
00223         double coef = frow->getEntry(rownum);
00224 
00225         frow->removeEntry(rownum);
00226 
00227         fei::CSVec csrow(*mrow);
00228 
00229         std::vector<int>& indices = csrow.indices();
00230         std::vector<double>& coefs = csrow.coefs();
00231 
00232         size_t rowlen = csrow.size();
00233         for(size_t ii=0; ii<rowlen; ++ii) {
00234           coefs[ii] *= coef;
00235         }
00236 
00237         frow->addEntries(rowlen, &coefs[0], &indices[0]);
00238         foundCoupling = true;
00239       }
00240     }
00241 
00242     if (foundCoupling == true) ++levelsOfCoupling;
00243     else finished = true;
00244   }
00245 
00246   if (levelsOfCoupling > 1) {
00247     FEI_CERR << "fei::removeCouplings WARNING, levelsOfCoupling="
00248       << levelsOfCoupling << " (Each level of coupling means that a slave DOF "
00249       << "depends on a master which is itself a slave.) Behavior is not well "
00250       << "understood for levelsOfCoupling greater than 1..."<<FEI_ENDL;
00251   }
00252 
00253   return levelsOfCoupling;
00254 }
00255 
00256 //----------------------------------------------------------------------------
00257 void global_union(MPI_Comm comm,
00258                   const fei::FillableMat& localMatrix,
00259                   fei::FillableMat& globalUnionMatrix)
00260 {
00261   globalUnionMatrix = localMatrix;
00262 
00263   std::vector<int> localintdata;
00264   std::vector<double> localdoubledata;
00265 
00266   int localProc = fei::localProc(comm);
00267   int numProcs = fei::numProcs(comm);
00268 
00269   if (numProcs < 2) {
00270     return;
00271   }
00272 
00273   //first pack the local matrix into a pair of std::vector objects
00274 
00275   pack_FillableMat(localMatrix, localintdata, localdoubledata);
00276 
00277   //next use Allgatherv to place every processor's packed arrays onto every
00278   //other processor.
00279 
00280   std::vector<int> recvintdatalengths;
00281   std::vector<int> recvintdata;
00282   int err = fei::Allgatherv(comm, localintdata, recvintdatalengths, recvintdata);
00283   if (err != 0) {
00284     throw std::runtime_error("fei::impl_utils::global_union: Allgatherv-int failed.");
00285   }
00286 
00287   std::vector<int> recvdoubledatalengths;
00288   std::vector<double> recvdoubledata;
00289   err = fei::Allgatherv(comm, localdoubledata, recvdoubledatalengths, recvdoubledata);
00290   if (err != 0) {
00291     throw std::runtime_error("fei::impl_utils::global_union: Allgatherv-double failed.");
00292   }
00293 
00294   if (recvintdatalengths.size() != recvdoubledatalengths.size()) {
00295     throw std::runtime_error("fei::impl_utils::global_union: inconsistent lengths from Allgatherv");
00296   }
00297 
00298   //finally unpack the received arrays into matrix objects and combine them
00299   //into the globalUnionMatrix object.
00300   bool clearMatOnEntry = false;
00301   bool overwriteEntries = true;
00302 
00303   int ioffset = 0;
00304   int doffset = 0;
00305   for(size_t p=0; p<recvintdatalengths.size(); ++p) {
00306     int intlen = recvintdatalengths[p];
00307     int doublelen = recvdoubledatalengths[p];
00308 
00309     if (intlen > 1 && (int)p != localProc) {
00310       std::vector<int> intdata(&recvintdata[ioffset], &recvintdata[ioffset]+intlen);
00311       std::vector<double> doubledata(&recvdoubledata[doffset], &recvdoubledata[doffset]+doublelen);
00312 
00313       unpack_FillableMat(intdata, doubledata, globalUnionMatrix,
00314                          clearMatOnEntry, overwriteEntries);
00315     }
00316 
00317     ioffset += intlen;
00318     doffset += doublelen;
00319   }
00320 }
00321 
00322 //----------------------------------------------------------------------------
00323 void global_union(MPI_Comm comm,
00324                   const fei::CSVec& localVec, fei::CSVec& globalUnionVec)
00325 {
00326   globalUnionVec.indices().clear();
00327   globalUnionVec.coefs().clear();
00328 
00329   const std::vector<int>& localindices = localVec.indices();
00330   const std::vector<double>& localcoefs = localVec.coefs();
00331 
00332   std::vector<int> localintdata;
00333   if (localindices.size() > 0) {
00334     localintdata.assign(&localindices[0], &localindices[0]+localindices.size());
00335   }
00336 
00337   std::vector<double> localdoubledata;
00338   if (localcoefs.size() > 0) {
00339     localdoubledata.assign(&localcoefs[0], &localcoefs[0]+localcoefs.size());
00340   }
00341 
00342   //use Allgatherv to place every processor's arrays onto every
00343   //other processor.
00344 
00345   std::vector<int> recvintdatalengths;
00346   std::vector<int> recvintdata;
00347   int err = fei::Allgatherv(comm, localintdata, recvintdatalengths, recvintdata);
00348   if (err != 0) {
00349     throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-int failed.");
00350   }
00351 
00352   std::vector<int> recvdoubledatalengths;
00353   std::vector<double> recvdoubledata;
00354   err = fei::Allgatherv(comm, localdoubledata, recvdoubledatalengths, recvdoubledata);
00355   if (err != 0) {
00356     throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-double failed.");
00357   }
00358 
00359   if (recvintdatalengths.size() != recvdoubledatalengths.size()) {
00360     throw std::runtime_error("snl_fei::globalUnion csvec: inconsistent lengths from Allgatherv");
00361   }
00362 
00363   //now unpack the received arrays into the globalUnionVec object.
00364 
00365   unsigned len = recvintdata.size();
00366   if (len > 0) {
00367     int* recvintPtr = &recvintdata[0];
00368     double* recvdoublePtr = &recvdoubledata[0];
00369 
00370     for(unsigned i=0; i<len; ++i) {
00371       fei::put_entry(globalUnionVec, recvintPtr[i], recvdoublePtr[i]);
00372     }
00373   }
00374 }
00375 
00376 //----------------------------------------------------------------------------
00377 void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSRMat& mat)
00378 {
00379   fei::SparseRowGraph& srg = mat.getGraph();
00380 
00381   std::vector<int>& rowNumbers = srg.rowNumbers;
00382   for(size_t i=0; i<rowNumbers.size(); ++i) {
00383     rowNumbers[i] = reducer.translateToReducedEqn(rowNumbers[i]);
00384   }
00385 
00386   std::vector<int>& colIndices = srg.packedColumnIndices;
00387   for(size_t i=0; i<colIndices.size(); ++i) {
00388     colIndices[i] = reducer.translateToReducedEqn(colIndices[i]);
00389   }
00390 }
00391 
00392 //----------------------------------------------------------------------------
00393 void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSVec& vec)
00394 {
00395   std::vector<int>& indices = vec.indices();
00396   for(size_t i=0; i<indices.size(); ++i) {
00397     indices[i] = reducer.translateToReducedEqn(indices[i]);
00398   }
00399 }
00400 
00401 //----------------------------------------------------------------------------
00402 void add_to_graph(const fei::CSRMat& inmat, fei::Graph& graph)
00403 {
00404   const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
00405   const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
00406   const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
00407 
00408   for(size_t i=0; i<rowNumbers.size(); ++i) {
00409     int row = rowNumbers[i];
00410     int offset = rowOffsets[i];
00411     int rowlen = rowOffsets[i+1]-offset;
00412     const int* indices = &pckColInds[offset];
00413 
00414     if (graph.addIndices(row, rowlen, indices) != 0) {
00415       throw std::runtime_error("fei::impl_utils::add_to_graph ERROR in graph.addIndices.");
00416     }
00417   }
00418 }
00419 
00420 //----------------------------------------------------------------------------
00421 void add_to_matrix(const fei::CSRMat& inmat, bool sum_into, fei::Matrix& matrix)
00422 {
00423   const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
00424   const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
00425   const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
00426   const std::vector<double>& pckCoefs = inmat.getPackedCoefs();
00427 
00428   for(size_t i=0; i<rowNumbers.size(); ++i) {
00429     int row = rowNumbers[i];
00430     int offset = rowOffsets[i];
00431     int rowlen = rowOffsets[i+1]-offset;
00432     const int* indices = &pckColInds[offset];
00433     const double* coefs = &pckCoefs[offset];
00434 
00435     if (sum_into) {
00436       if (matrix.sumIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
00437         throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
00438       }
00439     }
00440     else {
00441       if (matrix.copyIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
00442         throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.copyIn.");
00443       }
00444     }
00445   }
00446 }
00447 
00448 }//namespace impl_utils
00449 }//namespace fei
00450 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends
Generated on Wed Apr 13 10:08:23 2011 for FEI by  doxygen 1.6.3