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