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