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 CSVec* v = iter->second;
00068     const std::vector<int>& v_ind = v->indices();
00069     const std::vector<double>& v_coef = v->coefs();
00070     for(size_t i=0; i<v_ind.size(); ++i) {
00071       colind_ptr[offset] = v_ind[i];
00072       coef_ptr[offset++] = v_coef[i];
00073     }
00074   }
00075 
00076   return *this;
00077 }
00078 
00079 CSRMat&
00080 CSRMat::operator+=(const CSRMat& src)
00081 {
00082   FillableMat tmp;
00083   add_CSRMat_to_FillableMat(*this, tmp);
00084   add_CSRMat_to_FillableMat(src, tmp);
00085   *this = tmp;
00086   return *this;
00087 }
00088 
00089 bool
00090 CSRMat::operator==(const CSRMat& rhs) const
00091 {
00092   if (getGraph() != rhs.getGraph()) return false;
00093   return getPackedCoefs() == rhs.getPackedCoefs();
00094 }
00095 
00096 bool
00097 CSRMat::operator!=(const CSRMat& rhs) const
00098 {
00099   return !(*this == rhs);
00100 }
00101 
00102 void multiply_CSRMat_CSVec(const CSRMat& A, const CSVec& x, CSVec& y)
00103 {
00104   //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
00105 
00106   const std::vector<int>& rows = A.getGraph().rowNumbers;
00107   const int* rowoffs = &(A.getGraph().rowOffsets[0]);
00108   const std::vector<int>& colinds = A.getGraph().packedColumnIndices;
00109   const double* Acoef = &(A.getPackedCoefs()[0]);
00110 
00111   const std::vector<int>& xind = x.indices();
00112   const std::vector<double>& xcoef = x.coefs();
00113 
00114   const double* xcoef_ptr = &xcoef[0];
00115   const int* xind_ptr = &xind[0];
00116   int xlen = xcoef.size();
00117 
00118   std::vector<int>& yind = y.indices();
00119   std::vector<double>& ycoef = y.coefs();
00120 
00121   unsigned nrows = A.getNumRows();
00122 
00123   yind.resize(nrows);
00124   ycoef.resize(nrows);
00125 
00126   int* yind_ptr = &yind[0];
00127   double* ycoef_ptr = &ycoef[0];
00128 
00129   int jbeg = *rowoffs++;
00130   for(unsigned i=0; i<nrows; ++i) {
00131     int jend = *rowoffs++;
00132 
00133     double sum = 0.0;
00134     while(jbeg<jend) {
00135       int xoff = fei::binarySearch(colinds[jbeg], xind_ptr, xlen);
00136 
00137       if (xoff > -1) {
00138         sum += Acoef[jbeg]*xcoef_ptr[xoff];
00139       }
00140       ++jbeg;
00141     }
00142 
00143     yind_ptr[i] = rows[i];
00144     ycoef_ptr[i] = sum;
00145   }
00146 }
00147 
00148 void multiply_trans_CSRMat_CSVec(const CSRMat& A, const CSVec& x, CSVec& y)
00149 {
00150   const std::vector<int>& rows = A.getGraph().rowNumbers;
00151   const int* rowoffs = &(A.getGraph().rowOffsets[0]);
00152   const int* colinds = &(A.getGraph().packedColumnIndices[0]);
00153   const double* Acoef = &(A.getPackedCoefs()[0]);
00154 
00155   const std::vector<int>& xind = x.indices();
00156   const std::vector<double>& xcoef = x.coefs();
00157 
00158   const double* xcoef_ptr = &xcoef[0];
00159 
00160   unsigned nrows = A.getNumRows();
00161 
00162   std::vector<int> offsets;
00163   fei::impl_utils::find_offsets(rows, xind, offsets);
00164   const int* offsetsptr = &offsets[0];
00165 
00166   fei::CSVec fy;
00167 
00168   int jbeg = *rowoffs++;
00169   for(unsigned i=0; i<nrows; ++i) {
00170     int jend = *rowoffs++;
00171 
00172     int xoff = offsetsptr[i];
00173     if (xoff < 0) {
00174       jbeg = jend;
00175       continue;
00176     }
00177 
00178     double xcoeff = xcoef_ptr[xoff];
00179 
00180     while(jbeg<jend) {
00181       add_entry(fy, colinds[jbeg],Acoef[jbeg]*xcoeff);
00182       ++jbeg;
00183     }
00184   }
00185 
00186   y = fy;
00187 }
00188 
00189 void multiply_CSRMat_CSRMat(const CSRMat& A, const CSRMat& B, CSRMat& C,
00190                             bool storeResultZeros)
00191 {
00192   //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
00193 
00194   fei::FillableMat fc;
00195 
00196   const std::vector<int>& Arows = A.getGraph().rowNumbers;
00197   const std::vector<int>& Brows = B.getGraph().rowNumbers;
00198   if (Arows.size() < 1 || Brows.size() < 1) {
00199     C = fc;
00200     return;
00201   }
00202   const int* Arowoffs = &(A.getGraph().rowOffsets[0]);
00203   const int* Acols = &(A.getGraph().packedColumnIndices[0]);
00204   const double* Acoefs = &(A.getPackedCoefs()[0]);
00205 
00206   const int* Browoffs = &(B.getGraph().rowOffsets[0]);
00207   const std::vector<int>& Bcols = B.getGraph().packedColumnIndices;
00208   const double* Bcoefs = &(B.getPackedCoefs()[0]);
00209 
00210   static double fei_min = std::numeric_limits<double>::min();
00211 
00212   int jbeg = *Arowoffs++;
00213   for(size_t i=0; i<Arows.size(); ++i) {
00214     int row = Arows[i];
00215     int jend = *Arowoffs++;
00216 
00217     fei::CSVec* fc_row = NULL;
00218     if (storeResultZeros) {
00219       fc_row = fc.create_or_getRow(row);
00220     }
00221     else {
00222       fc_row = fc.hasRow(row) ? fc.create_or_getRow(row) : NULL;
00223     }
00224 
00225     while(jbeg<jend) {
00226       ++jbeg;
00227       int Acol = *Acols++;
00228       double Acoef = *Acoefs++;
00229 
00230       int Brow_offset = fei::binarySearch(Acol, &Brows[0], Brows.size());
00231 
00232       if (Brow_offset < 0) {
00233         continue;
00234       }
00235 
00236       if (!storeResultZeros) {
00237         if (std::abs(Acoef) < fei_min) {
00238           continue;
00239         }
00240       }
00241 
00242       const int* Brow_cols = &(Bcols[Browoffs[Brow_offset]]);
00243       const double* Brow_coefs = &(Bcoefs[Browoffs[Brow_offset]]);
00244       int Brow_len = Browoffs[Brow_offset+1]-Browoffs[Brow_offset];
00245 
00246       for(int k=0; k<Brow_len; ++k) {
00247         double resultCoef = Acoef*Brow_coefs[k];
00248         int resultCol = Brow_cols[k];
00249 
00250         if (!storeResultZeros) {
00251           if (std::abs(resultCoef) < fei_min) {
00252             continue;
00253           }
00254         }
00255 
00256         if (fc_row == NULL) {
00257           fc_row = fc.create_or_getRow(row);
00258         }
00259 
00260         add_entry(*fc_row, resultCol, resultCoef);
00261       }
00262     }
00263   }
00264 
00265   C = fc;
00266 }
00267 
00268 void multiply_trans_CSRMat_CSRMat(const CSRMat& A, const CSRMat& B, CSRMat& C,
00269                                   bool storeResultZeros)
00270 {
00271   //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
00272 
00273   fei::FillableMat fc;
00274 
00275   const std::vector<int>& Arows = A.getGraph().rowNumbers;
00276   const std::vector<int>& Brows = B.getGraph().rowNumbers;
00277   if (Arows.size() < 1 || Brows.size() < 1) {
00278     C = fc;
00279     return;
00280   }
00281 
00282   const size_t numArows = Arows.size();
00283   const int* Arowoffs = &(A.getGraph().rowOffsets[0]);
00284   const int* Acols = &(A.getGraph().packedColumnIndices[0]);
00285   const double* Acoefs = &(A.getPackedCoefs()[0]);
00286 
00287   const int* Browoffs = &(B.getGraph().rowOffsets[0]);
00288   const std::vector<int>& Bcols = B.getGraph().packedColumnIndices;
00289   const double* Bcoefs = &(B.getPackedCoefs()[0]);
00290 
00291   std::vector<double> row_coefs;
00292 
00293   static double fei_min = std::numeric_limits<double>::min();
00294 
00295   std::vector<int> offsets;
00296   fei::impl_utils::find_offsets(Arows, Brows, offsets);
00297 
00298   int jbeg = *Arowoffs++;
00299   for(size_t i=0; i<numArows; ++i) {
00300     int jend = *Arowoffs++;
00301 
00302     int Brow_offset = offsets[i];
00303     if (Brow_offset < 0) {
00304       jbeg = jend;
00305       continue;
00306     }
00307 
00308     const int* Brow_cols = &(Bcols[Browoffs[Brow_offset]]);
00309     const double* Brow_coefs = &(Bcoefs[Browoffs[Brow_offset]]);
00310     int Brow_len = Browoffs[Brow_offset+1]-Browoffs[Brow_offset];
00311 
00312     if ((int)row_coefs.size() < Brow_len) row_coefs.resize(Brow_len*2);
00313     double* row_coefs_ptr = &row_coefs[0];
00314 
00315     while(jbeg<jend) {
00316       int Acol = Acols[jbeg];
00317       double Acoef = Acoefs[jbeg++];
00318 
00319       if (std::abs(Acoef) < fei_min && !storeResultZeros) {
00320         continue;
00321       }
00322 
00323       for(int k=0; k<Brow_len; ++k) {
00324         row_coefs_ptr[k] = Acoef*Brow_coefs[k];
00325       }
00326 
00327       fc.sumInRow(Acol, Brow_cols, row_coefs_ptr, Brow_len);
00328     }
00329   }
00330 
00331   C = fc;
00332 }
00333 
00334 void add_CSRMat_to_FillableMat(const CSRMat& csrm, FillableMat& fm)
00335 {
00336   const std::vector<int>& rows = csrm.getGraph().rowNumbers;
00337   const int* rowoffs = &(csrm.getGraph().rowOffsets[0]);
00338   const std::vector<int>& cols = csrm.getGraph().packedColumnIndices;
00339   const double* coefs = &(csrm.getPackedCoefs()[0]);
00340 
00341   for(size_t i=0; i<rows.size(); ++i) {
00342     int row = rows[i];
00343 
00344     for(int j=rowoffs[i]; j<rowoffs[i+1]; ++j) {
00345       fm.sumInCoef(row, cols[j], coefs[j]);
00346     }
00347   }
00348 }
00349 
00350 }//namespace fei
00351 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends