FEI Version of the Day
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_FillableMat.hpp>
00013 #include <fei_CSRMat.hpp>
00014 #include <fei_CSVec.hpp>
00015 #include <fei_Graph.hpp>
00016 #include <fei_Matrix.hpp>
00017 #include <fei_Reducer.hpp>
00018 
00019 namespace fei {
00020 namespace impl_utils {
00021 
00022 //----------------------------------------------------------------------------
00023 void find_offsets(const std::vector<int>& sources,
00024                   const std::vector<int>& targets,
00025                   std::vector<int>& offsets)
00026 {
00027   offsets.assign(sources.size(), -1);
00028 
00029   size_t ssize = sources.size(), tsize = targets.size();
00030   size_t si = 0, ti = 0;
00031 
00032   const int* sourcesptr = &sources[0];
00033   const int* targetsptr = &targets[0];
00034 
00035   while(si<ssize && ti<tsize) {
00036     if (sourcesptr[si] == targetsptr[ti]) {
00037       offsets[si++] = ti++;
00038       continue;
00039     }
00040 
00041     while(sourcesptr[si] < targetsptr[ti] && si<ssize) {
00042       ++si;
00043     }
00044 
00045     while(sourcesptr[si] > targetsptr[ti] && ti<tsize) {
00046       ++ti;
00047     }
00048   }
00049 }
00050 
00051 //----------------------------------------------------------------------------
00052 size_t num_bytes_FillableMat(const fei::FillableMat& mat)
00053 {
00054   int nrows = mat.getNumRows();
00055   int nnz = fei::count_nnz(mat);
00056 
00057   int num_chars_int = (2 + nrows*2 + nnz)*sizeof(int);
00058   int num_chars_double = nnz*sizeof(double);
00059 
00060   return num_chars_int + num_chars_double;
00061 }
00062 
00063 //----------------------------------------------------------------------------
00064 void pack_FillableMat(const fei::FillableMat& mat, 
00065                       char* buffer)
00066 {
00067   int nrows = mat.getNumRows();
00068   int nnz = fei::count_nnz(mat);
00069 
00070   int num_chars_int = (2 + nrows*2 + nnz)*sizeof(int);
00071 
00072   int* intdata = reinterpret_cast<int*>(buffer);
00073   double* doubledata = reinterpret_cast<double*>(buffer+num_chars_int);
00074 
00075   int ioffset = 0;
00076   int doffset = 0;
00077 
00078   intdata[ioffset++] = nrows;
00079   intdata[ioffset++] = nnz;
00080 
00081   int ioffsetcols = 2+nrows*2;
00082 
00083   fei::FillableMat::const_iterator
00084     r_iter = mat.begin(),
00085     r_end = mat.end();
00086 
00087   for(; r_iter!=r_end; ++r_iter) {
00088     int rowNumber = r_iter->first;
00089     const fei::CSVec* row = r_iter->second;
00090 
00091     intdata[ioffset++] = rowNumber;
00092     const int rowlen = row->size();
00093     intdata[ioffset++] = rowlen;
00094 
00095     const std::vector<int>& rowindices = row->indices();
00096     const std::vector<double>& rowcoefs = row->coefs();
00097     for(int i=0; i<rowlen; ++i) {
00098       intdata[ioffsetcols++] = rowindices[i];
00099       doubledata[doffset++] = rowcoefs[i];
00100     }
00101   }
00102 }
00103 
00104 //----------------------------------------------------------------------------
00105 void unpack_FillableMat(const char* buffer_begin, const char* buffer_end,
00106                         fei::FillableMat& mat,
00107                         bool clear_mat_on_entry,
00108                         bool overwrite_entries)
00109 {
00110   if (clear_mat_on_entry) {
00111     mat.clear();
00112   }
00113 
00114   if (buffer_end == buffer_begin) {
00115     return;
00116   }
00117 
00118   const int* intdata = reinterpret_cast<const int*>(buffer_begin);
00119   int ioffset = 0;
00120   int nrows = intdata[ioffset++];
00121   int nnz = intdata[ioffset++];
00122 
00123   int ioffsetcols = 2+nrows*2;
00124 
00125   int num_chars_int = (2+nrows*2 + nnz)*sizeof(int);
00126   const double* doubledata = reinterpret_cast<const double*>(buffer_begin+num_chars_int);
00127 
00128   int doffset = 0;
00129 
00130   for(int i=0; i<nrows; ++i) {
00131     int row = intdata[ioffset++];
00132     int rowlen = intdata[ioffset++];
00133 
00134     for(int j=0; j<rowlen; ++j) {
00135       int col = intdata[ioffsetcols++];
00136       double coef = doubledata[doffset++];
00137 
00138       if (overwrite_entries) {
00139         mat.putCoef(row, col, coef);
00140       }
00141       else {
00142         mat.sumInCoef(row, col, coef);
00143       }
00144     }
00145   }
00146 
00147   if (doffset != nnz) {
00148     throw std::runtime_error("fei::impl_utils::unpack_FillableMat: failed, sizes don't agree.");
00149   }
00150 }
00151 
00152 //----------------------------------------------------------------------------
00153 bool unpack_CSRMat(const char* buffer_begin, const char* buffer_end, fei::CSRMat& mat)
00154 {
00155   bool all_zeros = true;
00156   if (buffer_end == buffer_begin) {
00157     return all_zeros;
00158   }
00159 
00160   const int* intdata = reinterpret_cast<const int*>(buffer_begin);
00161   int ioffset = 0;
00162   int nrows = intdata[ioffset++];
00163   int nnz = intdata[ioffset++];
00164 
00165   fei::SparseRowGraph& srg = mat.getGraph();
00166   srg.rowNumbers.resize(nrows);
00167   srg.rowOffsets.resize(nrows+1);
00168   srg.packedColumnIndices.resize(nnz);
00169   std::vector<double>& packed_coefs = mat.getPackedCoefs();
00170   packed_coefs.resize(nnz);
00171 
00172   int ioffsetcols = 2+nrows*2;
00173 
00174   int num_chars_int = (2+nrows*2 + nnz)*sizeof(int);
00175   const double* doubledata = reinterpret_cast<const double*>(buffer_begin+num_chars_int);
00176 
00177   int doffset = 0;
00178 
00179   for(int i=0; i<nrows; ++i) {
00180     int row = intdata[ioffset++];
00181     int rowlen = intdata[ioffset++];
00182 
00183     srg.rowNumbers[i] = row;
00184     srg.rowOffsets[i] = doffset;
00185 
00186     for(int j=0; j<rowlen; ++j) {
00187       int col = intdata[ioffsetcols++];
00188       double coef = doubledata[doffset];
00189       if (coef != 0.0) all_zeros = false;
00190       srg.packedColumnIndices[doffset] = col;
00191       packed_coefs[doffset++] = coef;
00192     }
00193   }
00194   srg.rowOffsets[nrows] = nnz;
00195   return all_zeros;
00196 }
00197 
00198 size_t num_bytes_indices_coefs(const std::vector<int>& indices,
00199                         const std::vector<double>& coefs)
00200 {
00201   int num = indices.size();
00202   int num_chars_int = (1+num)*sizeof(int);
00203   int num_chars = num_chars_int + num*sizeof(double);
00204   return num_chars;
00205 }
00206 
00207 void pack_indices_coefs(const std::vector<int>& indices,
00208                         const std::vector<double>& coefs,
00209                         std::vector<char>& buffer,
00210                         bool resize_buffer)
00211 {
00212   if (indices.size() != coefs.size()) {
00213     throw std::runtime_error("fei::impl_utils::pack_indices_coefs failed, sizes don't match.");
00214   }
00215 
00216   int num = indices.size();
00217   int num_chars_int = (1+num)*sizeof(int);
00218   int num_chars = num_chars_int + num*sizeof(double);
00219   if (resize_buffer) {
00220     buffer.resize(num_chars);
00221   }
00222 
00223   int* intdata = reinterpret_cast<int*>(&buffer[0]);
00224   double* doubledata = reinterpret_cast<double*>(&buffer[0]+num_chars_int);
00225 
00226   int ioffset = 0;
00227   int doffset = 0;
00228   intdata[ioffset++] = num;
00229   for(int i=0; i<num; ++i) {
00230     intdata[ioffset++] = indices[i];
00231     doubledata[doffset++] = coefs[i];
00232   }
00233 }
00234 
00235 void unpack_indices_coefs(const std::vector<char>& buffer,
00236                           std::vector<int>& indices,
00237                           std::vector<double>& coefs)
00238 {
00239   if (buffer.size() == 0) return;
00240 
00241   const int* intdata = reinterpret_cast<const int*>(&buffer[0]);
00242   int ioffset = 0;
00243   int num = intdata[ioffset++];
00244   int num_chars_int = (1+num)*sizeof(int);
00245   const double* doubledata = reinterpret_cast<const double*>(&buffer[0]+num_chars_int);
00246 
00247   indices.resize(num);
00248   coefs.resize(num);
00249 
00250   int doffset = 0;
00251   for(int i=0; i<num; ++i) {
00252     indices[i] = intdata[ioffset++];
00253     coefs[i] = doubledata[doffset++];
00254   }
00255 }
00256 
00257 //----------------------------------------------------------------------------
00258 void separate_BC_eqns(const fei::FillableMat& mat,
00259                     std::vector<int>& bcEqns,
00260                     std::vector<double>& bcVals)
00261 {
00262   fei::FillableMat::const_iterator
00263     m_iter = mat.begin(),
00264     m_end = mat.end();
00265 
00266   for(; m_iter != m_end; ++m_iter) {
00267     int eqn = m_iter->first;
00268     const fei::CSVec* row = m_iter->second;
00269 
00270     std::vector<int>::iterator
00271       iter = std::lower_bound(bcEqns.begin(), bcEqns.end(), eqn);
00272     if (iter == bcEqns.end() || *iter != eqn) {
00273       size_t offset = iter - bcEqns.begin();
00274       bcEqns.insert(iter, eqn);
00275       std::vector<double>::iterator viter = bcVals.begin();
00276       viter += offset;
00277       try {
00278         double val = get_entry(*row, eqn);
00279         bcVals.insert(viter, val);
00280       }
00281       catch(...) {
00282         fei::console_out() << "separate_BC_eqns ERROR, failed to find coef for eqn " << eqn;
00283         throw;
00284       }
00285     }
00286   }
00287 }
00288 
00289 //----------------------------------------------------------------------------
00290 void create_col_to_row_map(const fei::FillableMat& mat,
00291                            std::multimap<int,int>& crmap)
00292 {
00293   crmap.clear();
00294 
00295   if (mat.getNumRows() == 0) return;
00296 
00297   std::vector<int> rowNumbers;
00298   fei::get_row_numbers(mat, rowNumbers);
00299 
00300   fei::FillableMat::const_iterator
00301     m_iter = mat.begin(),
00302     m_end = mat.end();
00303 
00304   for(; m_iter != m_end; ++m_iter) {
00305     int rowNum = m_iter->first;
00306     const fei::CSVec* rowvec = m_iter->second;
00307 
00308     const std::vector<int>& rowindices = rowvec->indices();
00309 
00310     for(size_t i=0; i<rowindices.size(); ++i) {
00311       int colNum = rowindices[i];
00312 
00313       crmap.insert(std::make_pair(colNum, rowNum));
00314     }
00315   }
00316 }
00317 
00318 //----------------------------------------------------------------------------
00319 int remove_couplings(fei::FillableMat& mat)
00320 {
00321   int levelsOfCoupling = 0;
00322 
00323   std::vector<int> rowNumbers;
00324   fei::get_row_numbers(mat, rowNumbers);
00325 
00326   bool finished = false;
00327   while(!finished) {
00328     std::multimap<int,int> crmap;
00329     create_col_to_row_map(mat, crmap);
00330 
00331     typedef std::multimap<int,int>::iterator MM_Iter;
00332 
00333     fei::FillableMat::iterator
00334       m_iter = mat.begin(),
00335       m_end = mat.end();
00336 
00337     bool foundCoupling = false;
00338     for(; m_iter != m_end; ++m_iter) {
00339       int rownum = m_iter->first;
00340       fei::CSVec* mrow = m_iter->second;
00341 
00342       //now find which rows contain 'rownum' as a column-index:
00343       std::pair<MM_Iter,MM_Iter> mmi = crmap.equal_range(rownum);
00344 
00345       //loop over the rows which contain 'rownum' as a column-index:
00346       for(MM_Iter cri = mmi.first; cri != mmi.second; ++cri) {
00347         int cri_row = cri->second;
00348 
00349         fei::CSVec* frow = mat.create_or_getRow(cri_row);
00350 
00351         double coef = get_entry(*frow,rownum);
00352 
00353         remove_entry(*frow, rownum);
00354 
00355         std::vector<int>& indices = mrow->indices();
00356         std::vector<double>& coefs = mrow->coefs();
00357 
00358         size_t rowlen = mrow->size();
00359         for(size_t ii=0; ii<rowlen; ++ii) {
00360           coefs[ii] *= coef;
00361         }
00362 
00363         add_entries(*frow, rowlen, &indices[0], &coefs[0]);
00364         foundCoupling = true;
00365       }
00366     }
00367 
00368     if (foundCoupling == true) ++levelsOfCoupling;
00369     else finished = true;
00370   }
00371 
00372   if (levelsOfCoupling > 1) {
00373     fei::console_out() << "fei::removeCouplings WARNING, levelsOfCoupling="
00374       << levelsOfCoupling << " (Each level of coupling means that a slave DOF "
00375       << "depends on a master which is itself a slave.) Behavior is not well "
00376       << "understood for levelsOfCoupling greater than 1..."<<FEI_ENDL;
00377   }
00378 
00379   return levelsOfCoupling;
00380 }
00381 
00382 //----------------------------------------------------------------------------
00383 void global_union(MPI_Comm comm,
00384                   const fei::FillableMat& localMatrix,
00385                   fei::FillableMat& globalUnionMatrix)
00386 {
00387   globalUnionMatrix = localMatrix;
00388 
00389   int localProc = fei::localProc(comm);
00390   int numProcs = fei::numProcs(comm);
00391 
00392   if (numProcs < 2) {
00393     return;
00394   }
00395 
00396   //first pack the local matrix into a pair of std::vector objects
00397 
00398   size_t num_bytes = num_bytes_FillableMat(localMatrix);
00399   std::vector<char> localchardata(num_bytes);
00400 
00401   pack_FillableMat(localMatrix, &localchardata[0]);
00402 
00403   //next use Allgatherv to place every processor's packed arrays onto every
00404   //other processor.
00405 
00406   std::vector<int> recvdatalengths;
00407   std::vector<char> recvdata;
00408   int err = fei::Allgatherv(comm, localchardata, recvdatalengths, recvdata);
00409   if (err != 0) {
00410     throw std::runtime_error("fei::impl_utils::global_union: Allgatherv-int failed.");
00411   }
00412 
00413   //finally unpack the received arrays into matrix objects and combine them
00414   //into the globalUnionMatrix object.
00415   bool clearMatOnEntry = false;
00416   bool overwriteEntries = true;
00417 
00418   int ioffset = 0;
00419   for(size_t p=0; p<recvdatalengths.size(); ++p) {
00420     int len = recvdatalengths[p];
00421 
00422     if (len > 1 && (int)p != localProc) {
00423       unpack_FillableMat(&recvdata[ioffset], &recvdata[ioffset]+len,
00424                 globalUnionMatrix, clearMatOnEntry, overwriteEntries);
00425     }
00426 
00427     ioffset += len;
00428   }
00429 }
00430 
00431 //----------------------------------------------------------------------------
00432 void global_union(MPI_Comm comm,
00433                   const fei::CSVec& localVec, fei::CSVec& globalUnionVec)
00434 {
00435   globalUnionVec.indices().clear();
00436   globalUnionVec.coefs().clear();
00437 
00438   const std::vector<int>& localindices = localVec.indices();
00439   const std::vector<double>& localcoefs = localVec.coefs();
00440 
00441   std::vector<int> localintdata;
00442   if (localindices.size() > 0) {
00443     localintdata.assign(&localindices[0], &localindices[0]+localindices.size());
00444   }
00445 
00446   std::vector<double> localdoubledata;
00447   if (localcoefs.size() > 0) {
00448     localdoubledata.assign(&localcoefs[0], &localcoefs[0]+localcoefs.size());
00449   }
00450 
00451   //use Allgatherv to place every processor's arrays onto every
00452   //other processor.
00453 
00454   std::vector<int> recvintdatalengths;
00455   std::vector<int> recvintdata;
00456   int err = fei::Allgatherv(comm, localintdata, recvintdatalengths, recvintdata);
00457   if (err != 0) {
00458     throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-int failed.");
00459   }
00460 
00461   std::vector<int> recvdoubledatalengths;
00462   std::vector<double> recvdoubledata;
00463   err = fei::Allgatherv(comm, localdoubledata, recvdoubledatalengths, recvdoubledata);
00464   if (err != 0) {
00465     throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-double failed.");
00466   }
00467 
00468   if (recvintdatalengths.size() != recvdoubledatalengths.size()) {
00469     throw std::runtime_error("snl_fei::globalUnion csvec: inconsistent lengths from Allgatherv");
00470   }
00471 
00472   //now unpack the received arrays into the globalUnionVec object.
00473 
00474   unsigned len = recvintdata.size();
00475   if (len > 0) {
00476     int* recvintPtr = &recvintdata[0];
00477     double* recvdoublePtr = &recvdoubledata[0];
00478 
00479     for(unsigned i=0; i<len; ++i) {
00480       fei::put_entry(globalUnionVec, recvintPtr[i], recvdoublePtr[i]);
00481     }
00482   }
00483 }
00484 
00485 //----------------------------------------------------------------------------
00486 void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSRMat& mat)
00487 {
00488   fei::SparseRowGraph& srg = mat.getGraph();
00489 
00490   std::vector<int>& rowNumbers = srg.rowNumbers;
00491   for(size_t i=0; i<rowNumbers.size(); ++i) {
00492     rowNumbers[i] = reducer.translateToReducedEqn(rowNumbers[i]);
00493   }
00494 
00495   std::vector<int>& colIndices = srg.packedColumnIndices;
00496   for(size_t i=0; i<colIndices.size(); ++i) {
00497     colIndices[i] = reducer.translateToReducedEqn(colIndices[i]);
00498   }
00499 }
00500 
00501 //----------------------------------------------------------------------------
00502 void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSVec& vec)
00503 {
00504   std::vector<int>& indices = vec.indices();
00505   for(size_t i=0; i<indices.size(); ++i) {
00506     indices[i] = reducer.translateToReducedEqn(indices[i]);
00507   }
00508 }
00509 
00510 //----------------------------------------------------------------------------
00511 void add_to_graph(const fei::CSRMat& inmat, fei::Graph& graph)
00512 {
00513   const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
00514   const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
00515   const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
00516 
00517   for(size_t i=0; i<rowNumbers.size(); ++i) {
00518     int row = rowNumbers[i];
00519     int offset = rowOffsets[i];
00520     int rowlen = rowOffsets[i+1]-offset;
00521     const int* indices = &pckColInds[offset];
00522 
00523     if (graph.addIndices(row, rowlen, indices) != 0) {
00524       throw std::runtime_error("fei::impl_utils::add_to_graph ERROR in graph.addIndices.");
00525     }
00526   }
00527 }
00528 
00529 //----------------------------------------------------------------------------
00530 void add_to_matrix(const fei::CSRMat& inmat, bool sum_into, fei::Matrix& matrix)
00531 {
00532   const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
00533   const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
00534   const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
00535   const std::vector<double>& pckCoefs = inmat.getPackedCoefs();
00536 
00537   for(size_t i=0; i<rowNumbers.size(); ++i) {
00538     int row = rowNumbers[i];
00539     int offset = rowOffsets[i];
00540     int rowlen = rowOffsets[i+1]-offset;
00541     const int* indices = &pckColInds[offset];
00542     const double* coefs = &pckCoefs[offset];
00543 
00544     if (sum_into) {
00545       if (matrix.sumIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
00546         throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
00547       }
00548     }
00549     else {
00550       if (matrix.copyIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
00551         throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.copyIn.");
00552       }
00553     }
00554   }
00555 }
00556 
00557 }//namespace impl_utils
00558 }//namespace fei
00559 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends