FEI Version of the Day
fei_CSRMat.cpp
00001 
00002 /*--------------------------------------------------------------------*/
00003 /*    Copyright 2005 Sandia Corporation.                              */
00004 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00005 /*    non-exclusive license for use of this work by or on behalf      */
00006 /*    of the U.S. Government.  Export of this program may require     */
00007 /*    a license from the United States Government.                    */
00008 /*--------------------------------------------------------------------*/
00009 
00010 #include "fei_CSRMat.hpp"
00011 #include <fei_impl_utils.hpp>
00012 #include "fei_ArrayUtils.hpp"
00013 #include <limits>
00014 #include <cmath>
00015 
00016 namespace fei {
00017 
00018 CSRMat::CSRMat()
00019  : srg_(),
00020    packedcoefs_()
00021 {
00022 }
00023 
00024 CSRMat::CSRMat(const FillableMat& fmat)
00025  : srg_(),
00026    packedcoefs_()
00027 {
00028   *this = fmat;
00029 }
00030 
00031 CSRMat::~CSRMat()
00032 {
00033 }
00034 
00035 CSRMat&
00036 CSRMat::operator=(const fei::FillableMat& src)
00037 {
00038   FillableMat::const_iterator iter = src.begin(), iter_end = src.end();
00039 
00040   unsigned nrows = src.getNumRows();
00041 
00042   srg_.rowNumbers.resize(nrows);
00043   srg_.rowOffsets.resize(nrows+1);
00044 
00045   unsigned nnz = 0;
00046   unsigned i = 0;
00047   for(; iter != iter_end; ++iter, ++i) {
00048     srg_.rowNumbers[i] = iter->first;
00049     srg_.rowOffsets[i] = nnz;
00050     nnz += iter->second->size();
00051   }
00052 
00053   srg_.rowOffsets[nrows] = nnz;
00054 
00055   srg_.packedColumnIndices.resize(nnz);
00056   packedcoefs_.resize(nnz);
00057 
00058   int* colind_ptr = (srg_.packedColumnIndices.size()
00059     ? &(srg_.packedColumnIndices[0]) : 0);
00060   double* coef_ptr = (packedcoefs_.size()
00061     ? &(packedcoefs_[0]) : 0);
00062 
00063   iter = src.begin();
00064 
00065   unsigned offset = 0;
00066   for(; iter != iter_end; ++iter) {
00067     const FillableVec* v = iter->second;
00068     FillableVec::const_iterator viter = v->begin(), vend = v->end();
00069     for(; viter != vend; ++viter) {
00070       colind_ptr[offset] = viter->first;
00071       coef_ptr[offset++] = viter->second;
00072     }
00073   }
00074 
00075   return *this;
00076 }
00077 
00078 CSRMat&
00079 CSRMat::operator+=(const CSRMat& src)
00080 {
00081   FillableMat tmp;
00082   add_CSRMat_to_FillableMat(*this, tmp);
00083   add_CSRMat_to_FillableMat(src, tmp);
00084   *this = tmp;
00085   return *this;
00086 }
00087 
00088 bool
00089 CSRMat::operator==(const CSRMat& rhs) const
00090 {
00091   if (getGraph() != rhs.getGraph()) return false;
00092   return getPackedCoefs() == rhs.getPackedCoefs();
00093 }
00094 
00095 bool
00096 CSRMat::operator!=(const CSRMat& rhs) const
00097 {
00098   return !(*this == rhs);
00099 }
00100 
00101 void multiply_CSRMat_CSVec(const CSRMat& A, const CSVec& x, CSVec& y)
00102 {
00103   //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
00104 
00105   const std::vector<int>& rows = A.getGraph().rowNumbers;
00106   const int* rowoffs = &(A.getGraph().rowOffsets[0]);
00107   const std::vector<int>& colinds = A.getGraph().packedColumnIndices;
00108   const double* Acoef = &(A.getPackedCoefs()[0]);
00109 
00110   const std::vector<int>& xind = x.indices();
00111   const std::vector<double>& xcoef = x.coefs();
00112 
00113   const double* xcoef_ptr = &xcoef[0];
00114   const int* xind_ptr = &xind[0];
00115   int xlen = xcoef.size();
00116 
00117   std::vector<int>& yind = y.indices();
00118   std::vector<double>& ycoef = y.coefs();
00119 
00120   unsigned nrows = A.getNumRows();
00121 
00122   yind.resize(nrows);
00123   ycoef.resize(nrows);
00124 
00125   int* yind_ptr = &yind[0];
00126   double* ycoef_ptr = &ycoef[0];
00127 
00128   int jbeg = *rowoffs++;
00129   for(unsigned i=0; i<nrows; ++i) {
00130     int jend = *rowoffs++;
00131 
00132     double sum = 0.0;
00133     while(jbeg<jend) {
00134       int xoff = fei::binarySearch(colinds[jbeg], xind_ptr, xlen);
00135 
00136       if (xoff > -1) {
00137         sum += Acoef[jbeg]*xcoef_ptr[xoff];
00138       }
00139       ++jbeg;
00140     }
00141 
00142     yind_ptr[i] = rows[i];
00143     ycoef_ptr[i] = sum;
00144   }
00145 }
00146 
00147 void multiply_trans_CSRMat_CSVec(const CSRMat& A, const CSVec& x, CSVec& y)
00148 {
00149   const std::vector<int>& rows = A.getGraph().rowNumbers;
00150   const int* rowoffs = &(A.getGraph().rowOffsets[0]);
00151   const int* colinds = &(A.getGraph().packedColumnIndices[0]);
00152   const double* Acoef = &(A.getPackedCoefs()[0]);
00153 
00154   const std::vector<int>& xind = x.indices();
00155   const std::vector<double>& xcoef = x.coefs();
00156 
00157   const double* xcoef_ptr = &xcoef[0];
00158 
00159   unsigned nrows = A.getNumRows();
00160 
00161   std::vector<int> offsets;
00162   fei::impl_utils::find_offsets(rows, xind, offsets);
00163   const int* offsetsptr = &offsets[0];
00164 
00165   fei::FillableVec fy;
00166 
00167   int jbeg = *rowoffs++;
00168   for(unsigned i=0; i<nrows; ++i) {
00169     int jend = *rowoffs++;
00170 
00171     int xoff = offsetsptr[i];
00172     if (xoff < 0) {
00173       jbeg = jend;
00174       continue;
00175     }
00176 
00177     double xcoeff = xcoef_ptr[xoff];
00178 
00179     while(jbeg<jend) {
00180       fy.addEntry(colinds[jbeg],Acoef[jbeg]*xcoeff);
00181       ++jbeg;
00182     }
00183   }
00184 
00185   y = fy;
00186 }
00187 
00188 void multiply_CSRMat_CSRMat(const CSRMat& A, const CSRMat& B, CSRMat& C,
00189                             bool storeResultZeros)
00190 {
00191   //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
00192 
00193   fei::FillableMat fc;
00194 
00195   const std::vector<int>& Arows = A.getGraph().rowNumbers;
00196   const std::vector<int>& Brows = B.getGraph().rowNumbers;
00197   if (Arows.size() < 1 || Brows.size() < 1) {
00198     C = fc;
00199     return;
00200   }
00201   const int* Arowoffs = &(A.getGraph().rowOffsets[0]);
00202   const int* Acols = &(A.getGraph().packedColumnIndices[0]);
00203   const double* Acoefs = &(A.getPackedCoefs()[0]);
00204 
00205   const int* Browoffs = &(B.getGraph().rowOffsets[0]);
00206   const std::vector<int>& Bcols = B.getGraph().packedColumnIndices;
00207   const double* Bcoefs = &(B.getPackedCoefs()[0]);
00208 
00209   static double fei_min = std::numeric_limits<double>::min();
00210 
00211   int jbeg = *Arowoffs++;
00212   for(size_t i=0; i<Arows.size(); ++i) {
00213     int row = Arows[i];
00214     int jend = *Arowoffs++;
00215 
00216     fei::FillableVec* fc_row = NULL;
00217     if (storeResultZeros) {
00218       fc_row = fc.create_or_getRow(row);
00219     }
00220     else {
00221       fc_row = fc.hasRow(row) ? fc.create_or_getRow(row) : NULL;
00222     }
00223 
00224     while(jbeg<jend) {
00225       ++jbeg;
00226       int Acol = *Acols++;
00227       double Acoef = *Acoefs++;
00228 
00229       int Brow_offset = fei::binarySearch(Acol, &Brows[0], Brows.size());
00230 
00231       if (Brow_offset < 0) {
00232         continue;
00233       }
00234 
00235       if (!storeResultZeros) {
00236         if (std::abs(Acoef) < fei_min) {
00237           continue;
00238         }
00239       }
00240 
00241       const int* Brow_cols = &(Bcols[Browoffs[Brow_offset]]);
00242       const double* Brow_coefs = &(Bcoefs[Browoffs[Brow_offset]]);
00243       int Brow_len = Browoffs[Brow_offset+1]-Browoffs[Brow_offset];
00244 
00245       for(int k=0; k<Brow_len; ++k) {
00246         double resultCoef = Acoef*Brow_coefs[k];
00247         int resultCol = Brow_cols[k];
00248 
00249         if (!storeResultZeros) {
00250           if (std::abs(resultCoef) < fei_min) {
00251             continue;
00252           }
00253         }
00254 
00255         if (fc_row == NULL) {
00256           fc_row = fc.create_or_getRow(row);
00257         }
00258 
00259         fc_row->addEntry(resultCol, resultCoef);
00260       }
00261     }
00262   }
00263 
00264   C = fc;
00265 }
00266 
00267 void multiply_trans_CSRMat_CSRMat(const CSRMat& A, const CSRMat& B, CSRMat& C,
00268                                   bool storeResultZeros)
00269 {
00270   //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
00271 
00272   fei::FillableMat fc;
00273 
00274   const std::vector<int>& Arows = A.getGraph().rowNumbers;
00275   const std::vector<int>& Brows = B.getGraph().rowNumbers;
00276   if (Arows.size() < 1 || Brows.size() < 1) {
00277     C = fc;
00278     return;
00279   }
00280 
00281   const size_t numArows = Arows.size();
00282   const int* Arowoffs = &(A.getGraph().rowOffsets[0]);
00283   const int* Acols = &(A.getGraph().packedColumnIndices[0]);
00284   const double* Acoefs = &(A.getPackedCoefs()[0]);
00285 
00286   const int* Browoffs = &(B.getGraph().rowOffsets[0]);
00287   const std::vector<int>& Bcols = B.getGraph().packedColumnIndices;
00288   const double* Bcoefs = &(B.getPackedCoefs()[0]);
00289 
00290   std::vector<double> row_coefs;
00291 
00292   static double fei_min = std::numeric_limits<double>::min();
00293 
00294   std::vector<int> offsets;
00295   fei::impl_utils::find_offsets(Arows, Brows, offsets);
00296 
00297   int jbeg = *Arowoffs++;
00298   for(size_t i=0; i<numArows; ++i) {
00299     int jend = *Arowoffs++;
00300 
00301     int Brow_offset = offsets[i];
00302     if (Brow_offset < 0) {
00303       jbeg = jend;
00304       continue;
00305     }
00306 
00307     const int* Brow_cols = &(Bcols[Browoffs[Brow_offset]]);
00308     const double* Brow_coefs = &(Bcoefs[Browoffs[Brow_offset]]);
00309     int Brow_len = Browoffs[Brow_offset+1]-Browoffs[Brow_offset];
00310 
00311     if ((int)row_coefs.size() < Brow_len) row_coefs.resize(Brow_len*2);
00312     double* row_coefs_ptr = &row_coefs[0];
00313 
00314     while(jbeg<jend) {
00315       int Acol = Acols[jbeg];
00316       double Acoef = Acoefs[jbeg++];
00317 
00318       if (std::abs(Acoef) < fei_min && !storeResultZeros) {
00319         continue;
00320       }
00321 
00322       for(int k=0; k<Brow_len; ++k) {
00323         row_coefs_ptr[k] = Acoef*Brow_coefs[k];
00324       }
00325 
00326       fc.sumInRow(Acol, Brow_cols, row_coefs_ptr, Brow_len);
00327     }
00328   }
00329 
00330   C = fc;
00331 }
00332 
00333 void add_CSRMat_to_FillableMat(const CSRMat& csrm, FillableMat& fm)
00334 {
00335   const std::vector<int>& rows = csrm.getGraph().rowNumbers;
00336   const int* rowoffs = &(csrm.getGraph().rowOffsets[0]);
00337   const std::vector<int>& cols = csrm.getGraph().packedColumnIndices;
00338   const double* coefs = &(csrm.getPackedCoefs()[0]);
00339 
00340   for(size_t i=0; i<rows.size(); ++i) {
00341     int row = rows[i];
00342 
00343     for(int j=rowoffs[i]; j<rowoffs[i+1]; ++j) {
00344       fm.sumInCoef(row, cols[j], coefs[j]);
00345     }
00346   }
00347 }
00348 
00349 }//namespace fei
00350 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends