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_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 pack_FillableMat(const fei::FillableMat& mat, 
00091                       std::vector<char>& buffer)
00092 {
00093   int nrows = mat.getNumRows();
00094   int nnz = fei::count_nnz(mat);
00095 
00096   int num_chars_int = (2 + nrows*2 + nnz)*sizeof(int);
00097   int num_chars_double = nnz*sizeof(double);
00098   buffer.resize(num_chars_int + num_chars_double);
00099 
00100   int* intdata = reinterpret_cast<int*>(&buffer[0]);
00101   double* doubledata = reinterpret_cast<double*>(&buffer[0]+num_chars_int);
00102 
00103   int ioffset = 0;
00104   int doffset = 0;
00105 
00106   intdata[ioffset++] = nrows;
00107   intdata[ioffset++] = nnz;
00108 
00109   int ioffsetcols = 2+nrows*2;
00110 
00111   fei::FillableMat::const_iterator
00112     r_iter = mat.begin(),
00113     r_end = mat.end();
00114 
00115   for(; r_iter!=r_end; ++r_iter) {
00116     int rowNumber = r_iter->first;
00117     const fei::FillableVec* row = r_iter->second;
00118 
00119     intdata[ioffset++] = rowNumber;
00120     intdata[ioffset++] = row->size();
00121 
00122     fei::FillableVec::const_iterator
00123       iter = row->begin(),
00124       iend = row->end();
00125     for(; iter!=iend; ++iter) {
00126       intdata[ioffsetcols++] = iter->first;
00127       doubledata[doffset++] = iter->second;
00128     }
00129   }
00130 }
00131 
00132 //----------------------------------------------------------------------------
00133 void unpack_FillableMat(const std::vector<int>& intdata,
00134                         const std::vector<double>& doubledata,
00135                         fei::FillableMat& mat,
00136                         bool clear_mat_on_entry,
00137                         bool overwrite_entries)
00138 {
00139   if (clear_mat_on_entry) {
00140     mat.clear();
00141   }
00142 
00143   if (intdata.size() < 1) {
00144     return;
00145   }
00146 
00147   int ioffset = 0;
00148   int doffset = 0;
00149 
00150   int nrows = intdata[ioffset++];
00151 
00152   for(int i=0; i<nrows; ++i) {
00153     int row = intdata[ioffset++];
00154     int rowlen = intdata[ioffset++];
00155 
00156     for(int j=0; j<rowlen; ++j) {
00157       int col = intdata[ioffset++];
00158       double coef = doubledata[doffset++];
00159 
00160       if (overwrite_entries) {
00161         mat.putCoef(row, col, coef);
00162       }
00163       else {
00164         mat.sumInCoef(row, col, coef);
00165       }
00166     }
00167   }
00168 }
00169 
00170 //----------------------------------------------------------------------------
00171 void unpack_FillableMat(const std::vector<char>& buffer,
00172                         fei::FillableMat& mat,
00173                         bool clear_mat_on_entry,
00174                         bool overwrite_entries)
00175 {
00176   if (clear_mat_on_entry) {
00177     mat.clear();
00178   }
00179 
00180   if (buffer.size() < 1) {
00181     return;
00182   }
00183 
00184   const int* intdata = reinterpret_cast<const int*>(&buffer[0]);
00185   int ioffset = 0;
00186   int nrows = intdata[ioffset++];
00187   int nnz = intdata[ioffset++];
00188 
00189   int ioffsetcols = 2+nrows*2;
00190 
00191   int num_chars_int = (2+nrows*2 + nnz)*sizeof(int);
00192   const double* doubledata = reinterpret_cast<const double*>(&buffer[0]+num_chars_int);
00193 
00194   int doffset = 0;
00195 
00196   for(int i=0; i<nrows; ++i) {
00197     int row = intdata[ioffset++];
00198     int rowlen = intdata[ioffset++];
00199 
00200     for(int j=0; j<rowlen; ++j) {
00201       int col = intdata[ioffsetcols++];
00202       double coef = doubledata[doffset++];
00203 
00204       if (overwrite_entries) {
00205         mat.putCoef(row, col, coef);
00206       }
00207       else {
00208         mat.sumInCoef(row, col, coef);
00209       }
00210     }
00211   }
00212 
00213   if (doffset != nnz) {
00214     throw std::runtime_error("fei::impl_utils::unpack_FillableMat: failed, sizes don't agree.");
00215   }
00216 }
00217 
00218 //----------------------------------------------------------------------------
00219 void unpack_CSRMat(const std::vector<char>& buffer, fei::CSRMat& mat)
00220 {
00221   if (buffer.size() < 1) {
00222     return;
00223   }
00224 
00225   const int* intdata = reinterpret_cast<const int*>(&buffer[0]);
00226   int ioffset = 0;
00227   int nrows = intdata[ioffset++];
00228   int nnz = intdata[ioffset++];
00229 
00230   fei::SparseRowGraph& srg = mat.getGraph();
00231   srg.rowNumbers.resize(nrows);
00232   srg.rowOffsets.resize(nrows+1);
00233   srg.packedColumnIndices.resize(nnz);
00234   std::vector<double>& packed_coefs = mat.getPackedCoefs();
00235   packed_coefs.resize(nnz);
00236 
00237   int ioffsetcols = 2+nrows*2;
00238 
00239   int num_chars_int = (2+nrows*2 + nnz)*sizeof(int);
00240   const double* doubledata = reinterpret_cast<const double*>(&buffer[0]+num_chars_int);
00241 
00242   int doffset = 0;
00243 
00244   for(int i=0; i<nrows; ++i) {
00245     int row = intdata[ioffset++];
00246     int rowlen = intdata[ioffset++];
00247 
00248     srg.rowNumbers[i] = row;
00249     srg.rowOffsets[i] = doffset;
00250 
00251     for(int j=0; j<rowlen; ++j) {
00252       int col = intdata[ioffsetcols++];
00253       double coef = doubledata[doffset];
00254 
00255       srg.packedColumnIndices[doffset] = col;
00256       packed_coefs[doffset++] = coef;
00257     }
00258   }
00259   srg.rowOffsets[nrows] = nnz;
00260 }
00261 
00262 void pack_indices_coefs(const std::vector<int>& indices,
00263                         const std::vector<double>& coefs,
00264                         std::vector<char>& buffer)
00265 {
00266   if (indices.size() != coefs.size()) {
00267     throw std::runtime_error("fei::impl_utils::pack_indices_coefs failed, sizes don't match.");
00268   }
00269 
00270   int num = indices.size();
00271   int num_chars_int = (1+num)*sizeof(int);
00272   int num_chars = num_chars_int + num*sizeof(double);
00273   buffer.resize(num_chars);
00274 
00275   int* intdata = reinterpret_cast<int*>(&buffer[0]);
00276   double* doubledata = reinterpret_cast<double*>(&buffer[0]+num_chars_int);
00277 
00278   int ioffset = 0;
00279   int doffset = 0;
00280   intdata[ioffset++] = num;
00281   for(int i=0; i<num; ++i) {
00282     intdata[ioffset++] = indices[i];
00283     doubledata[doffset++] = coefs[i];
00284   }
00285 }
00286 
00287 void unpack_indices_coefs(const std::vector<char>& buffer,
00288                           std::vector<int>& indices,
00289                           std::vector<double>& coefs)
00290 {
00291   if (buffer.size() == 0) return;
00292 
00293   const int* intdata = reinterpret_cast<const int*>(&buffer[0]);
00294   int ioffset = 0;
00295   int num = intdata[ioffset++];
00296   int num_chars_int = (1+num)*sizeof(int);
00297   const double* doubledata = reinterpret_cast<const double*>(&buffer[0]+num_chars_int);
00298 
00299   indices.resize(num);
00300   coefs.resize(num);
00301 
00302   int doffset = 0;
00303   for(int i=0; i<num; ++i) {
00304     indices[i] = intdata[ioffset++];
00305     coefs[i] = doubledata[doffset++];
00306   }
00307 }
00308 
00309 //----------------------------------------------------------------------------
00310 void separate_BC_eqns(const fei::FillableMat& mat,
00311                     std::vector<int>& bcEqns,
00312                     std::vector<double>& bcVals)
00313 {
00314   fei::FillableMat::const_iterator
00315     m_iter = mat.begin(),
00316     m_end = mat.end();
00317 
00318   for(; m_iter != m_end; ++m_iter) {
00319     int eqn = m_iter->first;
00320     const fei::FillableVec* row = m_iter->second;
00321 
00322     std::vector<int>::iterator
00323       iter = std::lower_bound(bcEqns.begin(), bcEqns.end(), eqn);
00324     if (iter == bcEqns.end() || *iter != eqn) {
00325       size_t offset = iter - bcEqns.begin();
00326       bcEqns.insert(iter, eqn);
00327       std::vector<double>::iterator viter = bcVals.begin();
00328       viter += offset;
00329       try {
00330         double val = row->getEntry(eqn);
00331         bcVals.insert(viter, val);
00332       }
00333       catch(...) {
00334         fei::console_out() << "separate_BC_eqns ERROR, failed to find coef for eqn " << eqn;
00335         throw;
00336       }
00337     }
00338   }
00339 }
00340 
00341 //----------------------------------------------------------------------------
00342 void create_col_to_row_map(const fei::FillableMat& mat,
00343                            std::multimap<int,int>& crmap)
00344 {
00345   crmap.clear();
00346 
00347   if (mat.getNumRows() == 0) return;
00348 
00349   std::vector<int> rowNumbers;
00350   fei::get_row_numbers(mat, rowNumbers);
00351 
00352   fei::FillableMat::const_iterator
00353     m_iter = mat.begin(),
00354     m_end = mat.end();
00355 
00356   for(; m_iter != m_end; ++m_iter) {
00357     int rowNum = m_iter->first;
00358     const fei::FillableVec* rowvec = m_iter->second;
00359 
00360     fei::FillableVec::const_iterator
00361       r_iter = rowvec->begin(),
00362       r_end = rowvec->end();
00363 
00364     for(; r_iter != r_end; ++r_iter) {
00365       int colNum = r_iter->first;
00366 
00367       crmap.insert(std::make_pair(colNum, rowNum));
00368     }
00369   }
00370 }
00371 
00372 //----------------------------------------------------------------------------
00373 int remove_couplings(fei::FillableMat& mat)
00374 {
00375   int levelsOfCoupling = 0;
00376 
00377   std::vector<int> rowNumbers;
00378   fei::get_row_numbers(mat, rowNumbers);
00379 
00380   bool finished = false;
00381   while(!finished) {
00382     std::multimap<int,int> crmap;
00383     create_col_to_row_map(mat, crmap);
00384 
00385     typedef std::multimap<int,int>::iterator MM_Iter;
00386 
00387     fei::FillableMat::iterator
00388       m_iter = mat.begin(),
00389       m_end = mat.end();
00390 
00391     bool foundCoupling = false;
00392     for(; m_iter != m_end; ++m_iter) {
00393       int rownum = m_iter->first;
00394       fei::FillableVec* mrow = m_iter->second;
00395 
00396       //now find which rows contain 'rownum' as a column-index:
00397       std::pair<MM_Iter,MM_Iter> mmi = crmap.equal_range(rownum);
00398 
00399       //loop over the rows which contain 'rownum' as a column-index:
00400       for(MM_Iter cri = mmi.first; cri != mmi.second; ++cri) {
00401         int cri_row = cri->second;
00402 
00403         fei::FillableVec* frow = mat.create_or_getRow(cri_row);
00404 
00405         double coef = frow->getEntry(rownum);
00406 
00407         frow->removeEntry(rownum);
00408 
00409         fei::CSVec csrow(*mrow);
00410 
00411         std::vector<int>& indices = csrow.indices();
00412         std::vector<double>& coefs = csrow.coefs();
00413 
00414         size_t rowlen = csrow.size();
00415         for(size_t ii=0; ii<rowlen; ++ii) {
00416           coefs[ii] *= coef;
00417         }
00418 
00419         frow->addEntries(rowlen, &coefs[0], &indices[0]);
00420         foundCoupling = true;
00421       }
00422     }
00423 
00424     if (foundCoupling == true) ++levelsOfCoupling;
00425     else finished = true;
00426   }
00427 
00428   if (levelsOfCoupling > 1) {
00429     fei::console_out() << "fei::removeCouplings WARNING, levelsOfCoupling="
00430       << levelsOfCoupling << " (Each level of coupling means that a slave DOF "
00431       << "depends on a master which is itself a slave.) Behavior is not well "
00432       << "understood for levelsOfCoupling greater than 1..."<<FEI_ENDL;
00433   }
00434 
00435   return levelsOfCoupling;
00436 }
00437 
00438 //----------------------------------------------------------------------------
00439 void global_union(MPI_Comm comm,
00440                   const fei::FillableMat& localMatrix,
00441                   fei::FillableMat& globalUnionMatrix)
00442 {
00443   globalUnionMatrix = localMatrix;
00444 
00445   std::vector<int> localintdata;
00446   std::vector<double> localdoubledata;
00447 
00448   int localProc = fei::localProc(comm);
00449   int numProcs = fei::numProcs(comm);
00450 
00451   if (numProcs < 2) {
00452     return;
00453   }
00454 
00455   //first pack the local matrix into a pair of std::vector objects
00456 
00457   pack_FillableMat(localMatrix, localintdata, localdoubledata);
00458 
00459   //next use Allgatherv to place every processor's packed arrays onto every
00460   //other processor.
00461 
00462   std::vector<int> recvintdatalengths;
00463   std::vector<int> recvintdata;
00464   int err = fei::Allgatherv(comm, localintdata, recvintdatalengths, recvintdata);
00465   if (err != 0) {
00466     throw std::runtime_error("fei::impl_utils::global_union: Allgatherv-int failed.");
00467   }
00468 
00469   std::vector<int> recvdoubledatalengths;
00470   std::vector<double> recvdoubledata;
00471   err = fei::Allgatherv(comm, localdoubledata, recvdoubledatalengths, recvdoubledata);
00472   if (err != 0) {
00473     throw std::runtime_error("fei::impl_utils::global_union: Allgatherv-double failed.");
00474   }
00475 
00476   if (recvintdatalengths.size() != recvdoubledatalengths.size()) {
00477     throw std::runtime_error("fei::impl_utils::global_union: inconsistent lengths from Allgatherv");
00478   }
00479 
00480   //finally unpack the received arrays into matrix objects and combine them
00481   //into the globalUnionMatrix object.
00482   bool clearMatOnEntry = false;
00483   bool overwriteEntries = true;
00484 
00485   int ioffset = 0;
00486   int doffset = 0;
00487   for(size_t p=0; p<recvintdatalengths.size(); ++p) {
00488     int intlen = recvintdatalengths[p];
00489     int doublelen = recvdoubledatalengths[p];
00490 
00491     if (intlen > 1 && (int)p != localProc) {
00492       std::vector<int> intdata(&recvintdata[ioffset], &recvintdata[ioffset]+intlen);
00493       std::vector<double> doubledata(&recvdoubledata[doffset], &recvdoubledata[doffset]+doublelen);
00494 
00495       unpack_FillableMat(intdata, doubledata, globalUnionMatrix,
00496                          clearMatOnEntry, overwriteEntries);
00497     }
00498 
00499     ioffset += intlen;
00500     doffset += doublelen;
00501   }
00502 }
00503 
00504 //----------------------------------------------------------------------------
00505 void global_union(MPI_Comm comm,
00506                   const fei::CSVec& localVec, fei::CSVec& globalUnionVec)
00507 {
00508   globalUnionVec.indices().clear();
00509   globalUnionVec.coefs().clear();
00510 
00511   const std::vector<int>& localindices = localVec.indices();
00512   const std::vector<double>& localcoefs = localVec.coefs();
00513 
00514   std::vector<int> localintdata;
00515   if (localindices.size() > 0) {
00516     localintdata.assign(&localindices[0], &localindices[0]+localindices.size());
00517   }
00518 
00519   std::vector<double> localdoubledata;
00520   if (localcoefs.size() > 0) {
00521     localdoubledata.assign(&localcoefs[0], &localcoefs[0]+localcoefs.size());
00522   }
00523 
00524   //use Allgatherv to place every processor's arrays onto every
00525   //other processor.
00526 
00527   std::vector<int> recvintdatalengths;
00528   std::vector<int> recvintdata;
00529   int err = fei::Allgatherv(comm, localintdata, recvintdatalengths, recvintdata);
00530   if (err != 0) {
00531     throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-int failed.");
00532   }
00533 
00534   std::vector<int> recvdoubledatalengths;
00535   std::vector<double> recvdoubledata;
00536   err = fei::Allgatherv(comm, localdoubledata, recvdoubledatalengths, recvdoubledata);
00537   if (err != 0) {
00538     throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-double failed.");
00539   }
00540 
00541   if (recvintdatalengths.size() != recvdoubledatalengths.size()) {
00542     throw std::runtime_error("snl_fei::globalUnion csvec: inconsistent lengths from Allgatherv");
00543   }
00544 
00545   //now unpack the received arrays into the globalUnionVec object.
00546 
00547   unsigned len = recvintdata.size();
00548   if (len > 0) {
00549     int* recvintPtr = &recvintdata[0];
00550     double* recvdoublePtr = &recvdoubledata[0];
00551 
00552     for(unsigned i=0; i<len; ++i) {
00553       fei::put_entry(globalUnionVec, recvintPtr[i], recvdoublePtr[i]);
00554     }
00555   }
00556 }
00557 
00558 //----------------------------------------------------------------------------
00559 void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSRMat& mat)
00560 {
00561   fei::SparseRowGraph& srg = mat.getGraph();
00562 
00563   std::vector<int>& rowNumbers = srg.rowNumbers;
00564   for(size_t i=0; i<rowNumbers.size(); ++i) {
00565     rowNumbers[i] = reducer.translateToReducedEqn(rowNumbers[i]);
00566   }
00567 
00568   std::vector<int>& colIndices = srg.packedColumnIndices;
00569   for(size_t i=0; i<colIndices.size(); ++i) {
00570     colIndices[i] = reducer.translateToReducedEqn(colIndices[i]);
00571   }
00572 }
00573 
00574 //----------------------------------------------------------------------------
00575 void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSVec& vec)
00576 {
00577   std::vector<int>& indices = vec.indices();
00578   for(size_t i=0; i<indices.size(); ++i) {
00579     indices[i] = reducer.translateToReducedEqn(indices[i]);
00580   }
00581 }
00582 
00583 //----------------------------------------------------------------------------
00584 void add_to_graph(const fei::CSRMat& inmat, fei::Graph& graph)
00585 {
00586   const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
00587   const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
00588   const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
00589 
00590   for(size_t i=0; i<rowNumbers.size(); ++i) {
00591     int row = rowNumbers[i];
00592     int offset = rowOffsets[i];
00593     int rowlen = rowOffsets[i+1]-offset;
00594     const int* indices = &pckColInds[offset];
00595 
00596     if (graph.addIndices(row, rowlen, indices) != 0) {
00597       throw std::runtime_error("fei::impl_utils::add_to_graph ERROR in graph.addIndices.");
00598     }
00599   }
00600 }
00601 
00602 //----------------------------------------------------------------------------
00603 void add_to_matrix(const fei::CSRMat& inmat, bool sum_into, fei::Matrix& matrix)
00604 {
00605   const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
00606   const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
00607   const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
00608   const std::vector<double>& pckCoefs = inmat.getPackedCoefs();
00609 
00610   for(size_t i=0; i<rowNumbers.size(); ++i) {
00611     int row = rowNumbers[i];
00612     int offset = rowOffsets[i];
00613     int rowlen = rowOffsets[i+1]-offset;
00614     const int* indices = &pckColInds[offset];
00615     const double* coefs = &pckCoefs[offset];
00616 
00617     if (sum_into) {
00618       if (matrix.sumIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
00619         throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
00620       }
00621     }
00622     else {
00623       if (matrix.copyIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
00624         throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.copyIn.");
00625       }
00626     }
00627   }
00628 }
00629 
00630 }//namespace impl_utils
00631 }//namespace fei
00632 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends