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 bool unpack_CSRMat(const std::vector<char>& buffer, fei::CSRMat& mat)
00255 {
00256   bool all_zeros = true;
00257   if (buffer.size() < 1) {
00258     return all_zeros;
00259   }
00260 
00261   const int* intdata = reinterpret_cast<const int*>(&buffer[0]);
00262   int ioffset = 0;
00263   int nrows = intdata[ioffset++];
00264   int nnz = intdata[ioffset++];
00265 
00266   fei::SparseRowGraph& srg = mat.getGraph();
00267   srg.rowNumbers.resize(nrows);
00268   srg.rowOffsets.resize(nrows+1);
00269   srg.packedColumnIndices.resize(nnz);
00270   std::vector<double>& packed_coefs = mat.getPackedCoefs();
00271   packed_coefs.resize(nnz);
00272 
00273   int ioffsetcols = 2+nrows*2;
00274 
00275   int num_chars_int = (2+nrows*2 + nnz)*sizeof(int);
00276   const double* doubledata = reinterpret_cast<const double*>(&buffer[0]+num_chars_int);
00277 
00278   int doffset = 0;
00279 
00280   for(int i=0; i<nrows; ++i) {
00281     int row = intdata[ioffset++];
00282     int rowlen = intdata[ioffset++];
00283 
00284     srg.rowNumbers[i] = row;
00285     srg.rowOffsets[i] = doffset;
00286 
00287     for(int j=0; j<rowlen; ++j) {
00288       int col = intdata[ioffsetcols++];
00289       double coef = doubledata[doffset];
00290       if (coef != 0.0) all_zeros = false;
00291       srg.packedColumnIndices[doffset] = col;
00292       packed_coefs[doffset++] = coef;
00293     }
00294   }
00295   srg.rowOffsets[nrows] = nnz;
00296   return all_zeros;
00297 }
00298 
00299 void pack_indices_coefs(const std::vector<int>& indices,
00300                         const std::vector<double>& coefs,
00301                         std::vector<char>& buffer)
00302 {
00303   if (indices.size() != coefs.size()) {
00304     throw std::runtime_error("fei::impl_utils::pack_indices_coefs failed, sizes don't match.");
00305   }
00306 
00307   int num = indices.size();
00308   int num_chars_int = (1+num)*sizeof(int);
00309   int num_chars = num_chars_int + num*sizeof(double);
00310   buffer.resize(num_chars);
00311 
00312   int* intdata = reinterpret_cast<int*>(&buffer[0]);
00313   double* doubledata = reinterpret_cast<double*>(&buffer[0]+num_chars_int);
00314 
00315   int ioffset = 0;
00316   int doffset = 0;
00317   intdata[ioffset++] = num;
00318   for(int i=0; i<num; ++i) {
00319     intdata[ioffset++] = indices[i];
00320     doubledata[doffset++] = coefs[i];
00321   }
00322 }
00323 
00324 void unpack_indices_coefs(const std::vector<char>& buffer,
00325                           std::vector<int>& indices,
00326                           std::vector<double>& coefs)
00327 {
00328   if (buffer.size() == 0) return;
00329 
00330   const int* intdata = reinterpret_cast<const int*>(&buffer[0]);
00331   int ioffset = 0;
00332   int num = intdata[ioffset++];
00333   int num_chars_int = (1+num)*sizeof(int);
00334   const double* doubledata = reinterpret_cast<const double*>(&buffer[0]+num_chars_int);
00335 
00336   indices.resize(num);
00337   coefs.resize(num);
00338 
00339   int doffset = 0;
00340   for(int i=0; i<num; ++i) {
00341     indices[i] = intdata[ioffset++];
00342     coefs[i] = doubledata[doffset++];
00343   }
00344 }
00345 
00346 //----------------------------------------------------------------------------
00347 void separate_BC_eqns(const fei::FillableMat& mat,
00348                     std::vector<int>& bcEqns,
00349                     std::vector<double>& bcVals)
00350 {
00351   fei::FillableMat::const_iterator
00352     m_iter = mat.begin(),
00353     m_end = mat.end();
00354 
00355   for(; m_iter != m_end; ++m_iter) {
00356     int eqn = m_iter->first;
00357     const fei::FillableVec* row = m_iter->second;
00358 
00359     std::vector<int>::iterator
00360       iter = std::lower_bound(bcEqns.begin(), bcEqns.end(), eqn);
00361     if (iter == bcEqns.end() || *iter != eqn) {
00362       size_t offset = iter - bcEqns.begin();
00363       bcEqns.insert(iter, eqn);
00364       std::vector<double>::iterator viter = bcVals.begin();
00365       viter += offset;
00366       try {
00367         double val = row->getEntry(eqn);
00368         bcVals.insert(viter, val);
00369       }
00370       catch(...) {
00371         fei::console_out() << "separate_BC_eqns ERROR, failed to find coef for eqn " << eqn;
00372         throw;
00373       }
00374     }
00375   }
00376 }
00377 
00378 //----------------------------------------------------------------------------
00379 void create_col_to_row_map(const fei::FillableMat& mat,
00380                            std::multimap<int,int>& crmap)
00381 {
00382   crmap.clear();
00383 
00384   if (mat.getNumRows() == 0) return;
00385 
00386   std::vector<int> rowNumbers;
00387   fei::get_row_numbers(mat, rowNumbers);
00388 
00389   fei::FillableMat::const_iterator
00390     m_iter = mat.begin(),
00391     m_end = mat.end();
00392 
00393   for(; m_iter != m_end; ++m_iter) {
00394     int rowNum = m_iter->first;
00395     const fei::FillableVec* rowvec = m_iter->second;
00396 
00397     fei::FillableVec::const_iterator
00398       r_iter = rowvec->begin(),
00399       r_end = rowvec->end();
00400 
00401     for(; r_iter != r_end; ++r_iter) {
00402       int colNum = r_iter->first;
00403 
00404       crmap.insert(std::make_pair(colNum, rowNum));
00405     }
00406   }
00407 }
00408 
00409 //----------------------------------------------------------------------------
00410 int remove_couplings(fei::FillableMat& mat)
00411 {
00412   int levelsOfCoupling = 0;
00413 
00414   std::vector<int> rowNumbers;
00415   fei::get_row_numbers(mat, rowNumbers);
00416 
00417   bool finished = false;
00418   while(!finished) {
00419     std::multimap<int,int> crmap;
00420     create_col_to_row_map(mat, crmap);
00421 
00422     typedef std::multimap<int,int>::iterator MM_Iter;
00423 
00424     fei::FillableMat::iterator
00425       m_iter = mat.begin(),
00426       m_end = mat.end();
00427 
00428     bool foundCoupling = false;
00429     for(; m_iter != m_end; ++m_iter) {
00430       int rownum = m_iter->first;
00431       fei::FillableVec* mrow = m_iter->second;
00432 
00433       //now find which rows contain 'rownum' as a column-index:
00434       std::pair<MM_Iter,MM_Iter> mmi = crmap.equal_range(rownum);
00435 
00436       //loop over the rows which contain 'rownum' as a column-index:
00437       for(MM_Iter cri = mmi.first; cri != mmi.second; ++cri) {
00438         int cri_row = cri->second;
00439 
00440         fei::FillableVec* frow = mat.create_or_getRow(cri_row);
00441 
00442         double coef = frow->getEntry(rownum);
00443 
00444         frow->removeEntry(rownum);
00445 
00446         fei::CSVec csrow(*mrow);
00447 
00448         std::vector<int>& indices = csrow.indices();
00449         std::vector<double>& coefs = csrow.coefs();
00450 
00451         size_t rowlen = csrow.size();
00452         for(size_t ii=0; ii<rowlen; ++ii) {
00453           coefs[ii] *= coef;
00454         }
00455 
00456         frow->addEntries(rowlen, &coefs[0], &indices[0]);
00457         foundCoupling = true;
00458       }
00459     }
00460 
00461     if (foundCoupling == true) ++levelsOfCoupling;
00462     else finished = true;
00463   }
00464 
00465   if (levelsOfCoupling > 1) {
00466     fei::console_out() << "fei::removeCouplings WARNING, levelsOfCoupling="
00467       << levelsOfCoupling << " (Each level of coupling means that a slave DOF "
00468       << "depends on a master which is itself a slave.) Behavior is not well "
00469       << "understood for levelsOfCoupling greater than 1..."<<FEI_ENDL;
00470   }
00471 
00472   return levelsOfCoupling;
00473 }
00474 
00475 //----------------------------------------------------------------------------
00476 void global_union(MPI_Comm comm,
00477                   const fei::FillableMat& localMatrix,
00478                   fei::FillableMat& globalUnionMatrix)
00479 {
00480   globalUnionMatrix = localMatrix;
00481 
00482   std::vector<int> localintdata;
00483   std::vector<double> localdoubledata;
00484 
00485   int localProc = fei::localProc(comm);
00486   int numProcs = fei::numProcs(comm);
00487 
00488   if (numProcs < 2) {
00489     return;
00490   }
00491 
00492   //first pack the local matrix into a pair of std::vector objects
00493 
00494   pack_FillableMat(localMatrix, localintdata, localdoubledata);
00495 
00496   //next use Allgatherv to place every processor's packed arrays onto every
00497   //other processor.
00498 
00499   std::vector<int> recvintdatalengths;
00500   std::vector<int> recvintdata;
00501   int err = fei::Allgatherv(comm, localintdata, recvintdatalengths, recvintdata);
00502   if (err != 0) {
00503     throw std::runtime_error("fei::impl_utils::global_union: Allgatherv-int failed.");
00504   }
00505 
00506   std::vector<int> recvdoubledatalengths;
00507   std::vector<double> recvdoubledata;
00508   err = fei::Allgatherv(comm, localdoubledata, recvdoubledatalengths, recvdoubledata);
00509   if (err != 0) {
00510     throw std::runtime_error("fei::impl_utils::global_union: Allgatherv-double failed.");
00511   }
00512 
00513   if (recvintdatalengths.size() != recvdoubledatalengths.size()) {
00514     throw std::runtime_error("fei::impl_utils::global_union: inconsistent lengths from Allgatherv");
00515   }
00516 
00517   //finally unpack the received arrays into matrix objects and combine them
00518   //into the globalUnionMatrix object.
00519   bool clearMatOnEntry = false;
00520   bool overwriteEntries = true;
00521 
00522   int ioffset = 0;
00523   int doffset = 0;
00524   for(size_t p=0; p<recvintdatalengths.size(); ++p) {
00525     int intlen = recvintdatalengths[p];
00526     int doublelen = recvdoubledatalengths[p];
00527 
00528     if (intlen > 1 && (int)p != localProc) {
00529       std::vector<int> intdata(&recvintdata[ioffset], &recvintdata[ioffset]+intlen);
00530       std::vector<double> doubledata(&recvdoubledata[doffset], &recvdoubledata[doffset]+doublelen);
00531 
00532       unpack_FillableMat(intdata, doubledata, globalUnionMatrix,
00533                          clearMatOnEntry, overwriteEntries);
00534     }
00535 
00536     ioffset += intlen;
00537     doffset += doublelen;
00538   }
00539 }
00540 
00541 //----------------------------------------------------------------------------
00542 void global_union(MPI_Comm comm,
00543                   const fei::CSVec& localVec, fei::CSVec& globalUnionVec)
00544 {
00545   globalUnionVec.indices().clear();
00546   globalUnionVec.coefs().clear();
00547 
00548   const std::vector<int>& localindices = localVec.indices();
00549   const std::vector<double>& localcoefs = localVec.coefs();
00550 
00551   std::vector<int> localintdata;
00552   if (localindices.size() > 0) {
00553     localintdata.assign(&localindices[0], &localindices[0]+localindices.size());
00554   }
00555 
00556   std::vector<double> localdoubledata;
00557   if (localcoefs.size() > 0) {
00558     localdoubledata.assign(&localcoefs[0], &localcoefs[0]+localcoefs.size());
00559   }
00560 
00561   //use Allgatherv to place every processor's arrays onto every
00562   //other processor.
00563 
00564   std::vector<int> recvintdatalengths;
00565   std::vector<int> recvintdata;
00566   int err = fei::Allgatherv(comm, localintdata, recvintdatalengths, recvintdata);
00567   if (err != 0) {
00568     throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-int failed.");
00569   }
00570 
00571   std::vector<int> recvdoubledatalengths;
00572   std::vector<double> recvdoubledata;
00573   err = fei::Allgatherv(comm, localdoubledata, recvdoubledatalengths, recvdoubledata);
00574   if (err != 0) {
00575     throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-double failed.");
00576   }
00577 
00578   if (recvintdatalengths.size() != recvdoubledatalengths.size()) {
00579     throw std::runtime_error("snl_fei::globalUnion csvec: inconsistent lengths from Allgatherv");
00580   }
00581 
00582   //now unpack the received arrays into the globalUnionVec object.
00583 
00584   unsigned len = recvintdata.size();
00585   if (len > 0) {
00586     int* recvintPtr = &recvintdata[0];
00587     double* recvdoublePtr = &recvdoubledata[0];
00588 
00589     for(unsigned i=0; i<len; ++i) {
00590       fei::put_entry(globalUnionVec, recvintPtr[i], recvdoublePtr[i]);
00591     }
00592   }
00593 }
00594 
00595 //----------------------------------------------------------------------------
00596 void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSRMat& mat)
00597 {
00598   fei::SparseRowGraph& srg = mat.getGraph();
00599 
00600   std::vector<int>& rowNumbers = srg.rowNumbers;
00601   for(size_t i=0; i<rowNumbers.size(); ++i) {
00602     rowNumbers[i] = reducer.translateToReducedEqn(rowNumbers[i]);
00603   }
00604 
00605   std::vector<int>& colIndices = srg.packedColumnIndices;
00606   for(size_t i=0; i<colIndices.size(); ++i) {
00607     colIndices[i] = reducer.translateToReducedEqn(colIndices[i]);
00608   }
00609 }
00610 
00611 //----------------------------------------------------------------------------
00612 void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSVec& vec)
00613 {
00614   std::vector<int>& indices = vec.indices();
00615   for(size_t i=0; i<indices.size(); ++i) {
00616     indices[i] = reducer.translateToReducedEqn(indices[i]);
00617   }
00618 }
00619 
00620 //----------------------------------------------------------------------------
00621 void add_to_graph(const fei::CSRMat& inmat, fei::Graph& graph)
00622 {
00623   const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
00624   const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
00625   const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
00626 
00627   for(size_t i=0; i<rowNumbers.size(); ++i) {
00628     int row = rowNumbers[i];
00629     int offset = rowOffsets[i];
00630     int rowlen = rowOffsets[i+1]-offset;
00631     const int* indices = &pckColInds[offset];
00632 
00633     if (graph.addIndices(row, rowlen, indices) != 0) {
00634       throw std::runtime_error("fei::impl_utils::add_to_graph ERROR in graph.addIndices.");
00635     }
00636   }
00637 }
00638 
00639 //----------------------------------------------------------------------------
00640 void add_to_matrix(const fei::CSRMat& inmat, bool sum_into, fei::Matrix& matrix)
00641 {
00642   const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
00643   const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
00644   const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
00645   const std::vector<double>& pckCoefs = inmat.getPackedCoefs();
00646 
00647   for(size_t i=0; i<rowNumbers.size(); ++i) {
00648     int row = rowNumbers[i];
00649     int offset = rowOffsets[i];
00650     int rowlen = rowOffsets[i+1]-offset;
00651     const int* indices = &pckColInds[offset];
00652     const double* coefs = &pckCoefs[offset];
00653 
00654     if (sum_into) {
00655       if (matrix.sumIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
00656         throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
00657       }
00658     }
00659     else {
00660       if (matrix.copyIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
00661         throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.copyIn.");
00662       }
00663     }
00664   }
00665 }
00666 
00667 }//namespace impl_utils
00668 }//namespace fei
00669 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends