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