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