largeMatrix/denseStorage/RowDenseStorage.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
  \file RowDenseStorage.cpp
  \authors D. Martin, E. Lunéville, NGUYEN Manh Ha, N. Kielbasiewicz
  \since 21 jun 2011
  \date 21 feb 2013

  \brief Implementation of xlifepp::RowDenseStorage class and functionnalities
*/

#include "RowDenseStorage.hpp"
#include "utils.h"
#include <functional>

namespace xlifepp
{
/*---------------------------------------------------------------------------------------
  Constructors, destructor
  ---------------------------------------------------------------------------------------*/
// default constructor
RowDenseStorage::RowDenseStorage(string_t id)
  : DenseStorage(_row, id) {}

// constructor for square structure (number of rows)
RowDenseStorage::RowDenseStorage(number_t nr, string_t id)
  : DenseStorage(_row, nr, id) {}

// constructor for rectangular structure (number of rows/columns)
RowDenseStorage::RowDenseStorage(number_t nr, number_t nc, string_t id)
  : DenseStorage(_row, nr, nc, id) {}

// // constructor from a SymDenseStorage
// RowDenseStorage::RowDenseStorage(const SymDenseStorage& sd)
// : DenseStorage(_row, sd.numberOfRows()) {}

// // constructor from a SymmDenseStorage
// RowDenseStorage::RowDenseStorage(const DualDenseStorage& dd)
// : DenseStorage(_row, dd.numberOfRows(), dd.numberOfColumns()) {}

/*--------------------------------------------------------------------------------
 access operator to position of (i,j) coefficients (return 0 if outside)
 This DIRECT access to (i,j) pointer SHOULD NOT be used in matrix computation
 where SEQUENTIAL access MUST ALWAYS be used and is only provided for tests.
 SymType is not used !
--------------------------------------------------------------------------------*/
number_t RowDenseStorage::pos(number_t i, number_t j, SymType s) const
{
  if(i == 0 || i > nbRows_ || j == 0 || j > nbCols_) { return 0; }
  else { return ((i - 1) * nbCols_ + j); }
}

// access to submatrix adresses, useful for assembling matrices
// return the vector of adresses (relative number in storage) stored by row
void RowDenseStorage::positions(const std::vector<number_t>& rows, const std::vector<number_t>& cols,
                                std::vector<number_t>& pos, bool errorOn, SymType sym) const
{
  number_t nbp = rows.size() * cols.size();
  if(pos.size() != nbp) { pos.resize(nbp); }
  std::vector<number_t>::iterator it_p = pos.begin();
  for(std::vector<number_t>::const_iterator it_r = rows.begin(); it_r != rows.end(); it_r++)
    for(std::vector<number_t>::const_iterator it_c = cols.begin(); it_c != cols.end(); it_c++, it_p++)
    { *it_p = nbCols_ * (*it_r - 1) + *it_c; }
}

// get (row indices, adress) of col c in set [r1,r2] ## SymType not used ##
std::vector<std::pair<number_t, number_t> > RowDenseStorage::getCol(SymType s, number_t c, number_t r1, number_t r2) const
{
    number_t nbr=r2;
    if(nbr==0) nbr=nbRows_;
    std::vector<std::pair<number_t, number_t> > rows(nbr-r1+1);
    std::vector<std::pair<number_t, number_t> >::iterator itrow=rows.begin();
    for(number_t i=r1;i<=nbr;i++,itrow++) *itrow=std::make_pair(i,(i - 1) * nbCols_+c);
    return rows;
}

// get (col indices, adress) of row r in set [c1,c2] ## SymType not used ##
std::vector<std::pair<number_t, number_t> > RowDenseStorage::getRow(SymType s, number_t r, number_t c1, number_t c2) const
{
    number_t nbc=c2;
    if(nbc==0) nbc=nbCols_;
    std::vector<std::pair<number_t, number_t> > cols(nbc-c1+1);
    std::vector<std::pair<number_t, number_t> >::iterator itcol=cols.begin();
    number_t a=(r - 1) * nbCols_+c1;
    for(number_t j=c1;j<=nbc;j++, a++, itcol++) *itcol=std::make_pair(j,a);
    return cols;
}

// create a new scalar RowDense storage from current RowDense storage and submatrix sizes
RowDenseStorage* RowDenseStorage::toScalar(dimen_t nbr, dimen_t nbc)
{
    return new RowDenseStorage(nbr*nbRows_, nbc*nbCols_, stringId+"_scalar");
}

/*--------------------------------------------------------------------------------
 template Matrix x Vector & Vector x Matrix multiplications and specializations
 result vector has to be sized before
--------------------------------------------------------------------------------*/
template<typename M, typename V, typename R>
void RowDenseStorage::multMatrixVector(const std::vector<M>& m, V* vp, R* rp) const
{
  typename std::vector<M>::const_iterator itm = m.begin() + 1;
  V* vpe=vp+nbCols_;
  R* rpe=rp+nbRows_;
  rowMatrixVector(itm, vp, vpe, rp, rpe);
}

template<typename M, typename V, typename R>
void RowDenseStorage::multMatrixVector(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
  //trace_p->push("RowDenseStorage::multMatrixVector");
  typename std::vector<M>::const_iterator itm = m.begin() + 1;
  typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end();
  typename std::vector<R>::iterator it_rb = r.begin(), it_re = r.end();
  #ifndef XLIFEPP_WITH_OMP
  rowMatrixVector(itm, itvb, itve, it_rb, it_re);
  #else
  if(Environment::parallelOn()) parallelRowMatrixVector(itm, v, r);
  else rowMatrixVector(itm, itvb, itve, it_rb, it_re);
  #endif // XLIFEPP_WITH_OMP
  //trace_p->pop();
}

template<typename M, typename V, typename R>
void RowDenseStorage::multVectorMatrix(const std::vector<M>& m, V* vp, R* rp) const
{
  typename std::vector<M>::const_iterator itm = m.begin() + 1;
  V* vpe=vp+nbCols_;
  R* rpe=rp+nbRows_;
  rowVectorMatrix(itm, vp, vpe, rp, rpe);
}

template<typename M, typename V, typename R>
void RowDenseStorage::multVectorMatrix(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
  //trace_p->push("RowDenseStorage::multMatrixVector");
  typename std::vector<M>::const_iterator itm = m.begin() + 1;
  typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end();
  typename std::vector<R>::iterator it_rb = r.begin(), it_re = r.end();
  #ifndef XLIFEPP_WITH_OMP
  rowVectorMatrix(itm, itvb, itve, it_rb, it_re);
  #else
  if(Environment::parallelOn()) parallelMultVectorMatrix(m,v,r);
  else rowVectorMatrix(itm, itvb, itve, it_rb, it_re);
  #endif // XLIFEPP_WITH_OMP
  //trace_p->pop();
}

// parallel Matrix * Vector (specific)
template<typename M, typename V, typename R>
void RowDenseStorage::parallelMultVectorMatrix(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
     #ifdef XLIFEPP_WITH_OMP
     //trace_p->push("RowDenseStorage::parallelMultMatrixVector");
     number_t nr=nbOfRows(), nc=nbOfColumns();
     number_t nt = 1;
     #pragma omp parallel for lastprivate(nt)
      for (number_t i = 0; i < 1; i++)
        {nt = omp_get_num_threads();}

     number_t dr = nr/nt;
     if(dr==0 || nt==1 || !Environment::parallelOn()) //do not use parallel computation
     {
         typename std::vector<M>::const_iterator itm = m.begin() + 1;
         typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end();
         typename std::vector<R>::iterator it_rb = r.begin(), it_re = r.end();
         rowVectorMatrix(itm, itvb, itve, it_rb, it_re);
         //trace_p->pop();
         return;
     }
     std::vector<std::vector<R> > res(nt, std::vector<R>(nc, 0.*(v[0]*m[1])));  //nt temporary res vector
     #pragma omp parallel for
     for(number_t t=0; t<nt; ++t)  //partial product in parallel
     {
         typename std::vector<M>::const_iterator itm = m.begin() + ( 1 + t*dr*nc);
         typename std::vector<R>::iterator itrb = res[t].begin(), itre = res[t].end(), itr;
         typename std::vector<V>::const_iterator itvb = v.begin()+ t*dr, itve = itvb + dr, itv;
         if(t==nt-1) itve=v.end();
         for(itv = itvb; itv != itve; ++itv)
           for(itr = itrb; itr != itre; ++itr, ++itm)  *itr += *itm * *itv;
     }
     typename std::vector<R>::iterator itrb = r.begin(), itre = r.end(), itr;
     typename std::vector<R>::iterator itrp;
     for(itr = itrb; itr != itre; ++itr) { *itr *= 0; }
     for(number_t t=0; t<nt; ++t) //accumulate in serial
     {
         itrp = res[t].begin();
         for(itr=itrb; itr!=itre; ++itr, ++itrp) *itr+=*itrp;
     }
     //trace_p->pop();
     #endif // XLIFEPP_WITH_OMP
}

/*-----------------------------------------------------------------------
              triangular part matrix * vector
--------------------------------------------------------------------------*/

template<typename M, typename V, typename R>
void RowDenseStorage::lowerMatrixVector(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
    r.assign(nbRows_,R(0));  //to reset to 0
    typename std::vector<M>::const_iterator itmb=m.begin()+1, itm;
    typename std::vector<V>::const_iterator itv=v.begin();
    #ifdef XLIFEPP_WITH_OMP
    #pragma omp parallel for firstprivate(itv,itm) schedule(dynamic)
    #endif // XLIFEPP_WITH_OMP
    for(number_t i=0;i<nbRows_; i++)
       {
           itv=v.begin();
           itm=itmb+i*nbCols_;
           R ri=R(0);
           for(number_t j=0;j<=std::min(i,nbCols_-1);j++, ++itv, ++itm) ri+= *itm * *itv;
           r[i]=ri;
       }
}

template<typename M, typename V, typename R>
void RowDenseStorage::lowerD1MatrixVector(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
    r.assign(nbRows_,R(0));  //to reset to 0
    typename std::vector<M>::const_iterator itmb=m.begin()+1, itm;
    typename std::vector<V>::const_iterator itv=v.begin();
    #ifdef XLIFEPP_WITH_OMP
    #pragma omp parallel for firstprivate(itv,itm) schedule(dynamic)
    #endif // XLIFEPP_WITH_OMP
    for(number_t i=0;i<nbRows_; i++)
       {

           itv=v.begin();
           itm=itmb+i*nbCols_;
           R ri=R(0);
           for(number_t j=0;j<std::min(i,nbCols_);j++, ++itv, ++itm) ri+= *itm * *itv;
           if(i<nbCols_) ri+=*itv;
           r[i]=ri;
       }
}

template<typename M, typename V, typename R>
void RowDenseStorage::upperMatrixVector(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
    r.assign(nbRows_,R(0));  //to reset to 0
    typename std::vector<M>::const_iterator itmb=m.begin()+1, itm;
    typename std::vector<V>::const_iterator itv=v.begin();
    #ifdef XLIFEPP_WITH_OMP
    #pragma omp parallel for firstprivate(itv,itm) schedule(dynamic)
    #endif // XLIFEPP_WITH_OMP
    for(number_t i=0;i<std::min(nbRows_,nbCols_); i++)
       {
           itv=v.begin()+i;
           itm=itmb+(i*nbCols_+i);
           R ri=R(0);
           for(number_t j=i;j<nbCols_;j++, ++itv, ++itm) ri+= *itm * *itv;
           r[i]=ri;
       }
}

template<typename M, typename V, typename R>
void RowDenseStorage::upperD1MatrixVector(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
    r.assign(nbRows_,R(0));  //to reset to 0
    typename std::vector<M>::const_iterator itmb=m.begin()+1, itm;
    typename std::vector<V>::const_iterator itv=v.begin();
    #ifdef XLIFEPP_WITH_OMP
    #pragma omp parallel for firstprivate(itv,itm) schedule(dynamic)
    #endif // XLIFEPP_WITH_OMP
    for(number_t i=0;i<std::min(nbRows_,nbCols_); i++)
       {
           itv=v.begin()+i;
           R ri=*itv; itv++;
           itm=itmb+(i*nbCols_+i+1);
           for(number_t j=i+1;j<nbCols_;j++, ++itv, ++itm) ri+= *itm * *itv;
           r[i]=ri;
       }
}

template<typename M, typename V, typename R>
void RowDenseStorage::diagonalMatrixVector(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{

    if(nbRows_>nbCols_) r.assign(nbRows_,R(0));  //to reset to 0
    else r.resize(nbRows_);
    typename std::vector<M>::const_iterator itm=m.begin()+1;
    typename std::vector<V>::const_iterator itv=v.begin();
    typename std::vector<R>::iterator itr=r.begin();
    for(number_t i=0;i<std::min(nbRows_,nbCols_);i++, ++itv,++itr, itm+=nbCols_) *itr = *itm * *itv;
}


/*--------------------------------------------------------------------------------
 print RowDense matrix to ostream, according to type of values of matrix
 printScalarEntries is defined in DenseStorage class
--------------------------------------------------------------------------------*/
void RowDenseStorage::printEntries(std::ostream& os, const std::vector<real_t>& m, number_t vb, const SymType sym) const
{
  std::vector<real_t>::const_iterator itm = m.begin() + 1;
  printScalarEntries(itm, nbRows_, nbCols_, entriesPerRow,
                       entryWidth, entryPrec, "row", vb, os);
}

void RowDenseStorage::printEntries(std::ostream& os, const std::vector<complex_t>& m, number_t vb, const SymType sym) const
{
  std::vector<complex_t>::const_iterator itm = m.begin() + 1;
  printScalarEntries(itm, nbRows_, nbCols_, entriesPerRow / 2,
                     2 * entryWidth + 1, entryPrec, "row", vb, os);
}

void RowDenseStorage::printEntries(std::ostream& os, const std::vector< Vector<real_t> >& m, number_t vb, const SymType sym) const
{error("storage_novector", "ColDenseStorage::printEntries");}

void RowDenseStorage::printEntries(std::ostream& os, const std::vector< Vector<complex_t> >& m, number_t vb, const SymType sym) const
{error("storage_novector", "ColDenseStorage::printEntries");}

void RowDenseStorage::printEntries(std::ostream& os, const std::vector< Matrix<real_t> >& m, number_t vb, const SymType sym) const
{
  std::vector< Matrix<real_t> >::const_iterator itm = m.begin() + 1;
  printMatrixEntries(itm, nbRows_, nbCols_, "row", vb, os);
}
void RowDenseStorage::printEntries(std::ostream& os, const std::vector< Matrix<complex_t> >& m, number_t vb, const SymType sym) const
{
  std::vector< Matrix<complex_t> >::const_iterator itm = m.begin() + 1;
  printMatrixEntries(itm, nbRows_, nbCols_, "row", vb, os);
}

// Set the Diagonal with a specific value of DualDense and SymDense
template<typename T>
void RowDenseStorage::setDiagValueRowDense(std::vector<T>& v, const T m)
{
  typename std::vector<T>::iterator it = v.begin() + 1;
  number_t diagSize = diagonalSize() + 1;
  *it = m;
  for (number_t k = 1; k < (diagSize - 1); k++) { it += diagSize; (*it) = m;}
}

/*!
  Add two matrices
  \param m1 vector values_ of first matrix
  \param m2 vector values_ of second matrix
  \param r vector values_ of result matrix
 */
template<typename M1, typename M2, typename R>
void RowDenseStorage::addMatrixMatrix(const std::vector<M1>& m1, const std::vector<M2>& m2, std::vector<R>& r) const
{
  trace_p->push("RowDenseStorage::addMatrixMatrix");
  typename std::vector<M1>::const_iterator itm1 = m1.begin() + 1;
  typename std::vector<M2>::const_iterator itm2 = m2.begin() + 1;
  typename std::vector<R>::iterator it_rb = r.begin() + 1, it_re = r.end();
  sumMatrixMatrix(itm1, itm2, it_rb, it_re);
  trace_p->pop();
}

/*!
   Extract and convert matrix storage to UMFPack format (Matlab sparse matrix)
   \param m1 vector values_ current matrix
   \param colPtUmf vector column Pointer of UMFPack format
   \param rowIdxUmf vector row Index of UMFPack format
   \param resultUmf vector values of UMFPack format
 */
template<typename M1, typename Idx>
void RowDenseStorage::toUmfPack(const std::vector<M1>& m1, std::vector<Idx>& colPtUmf, std::vector<Idx>& rowIdxUmf, std::vector<M1>& resultUmf) const
{
    // Resize column pointer as dualCs's column pointer size
    colPtUmf.clear();
    colPtUmf.resize(nbCols_+1, Idx(0));

    // Because colPointer always begins with 0, so need to to count from the second position
    typename std::vector<Idx>::iterator itbColPtUmf = colPtUmf.begin(), itColPtUmf;

    typename std::vector<M1>::const_iterator itbval = m1.begin()+1;
    typename std::vector<M1>::const_iterator itval, iterange, itbrange;

    typename std::list<M1>::iterator itLVal;
    std::list<number_t>::iterator itLRow;

    std::list<M1> listValues;
    std::list<number_t> listRowIdx;

    std::vector<Idx> nbNzCol(nbCols_, Idx(0));
    typename std::vector<Idx>::const_iterator itbNzCol = nbNzCol.begin(), iteNzCol = nbNzCol.end(), itNzCol;

    for (int_t rowIdx = nbRows_-1; rowIdx  >= 0; --rowIdx) {
        itval = itbrange = itbval + rowIdx*nbCols_;
        iterange = itbrange + nbCols_;

        while (itval != iterange) {
            itval = std::find_if(itval, iterange, std::bind2nd(std::not_equal_to<M1>(), M1(0)));
            if (itval != iterange){
                number_t pos = itval - itbrange;

                ++nbNzCol[pos];
                itNzCol = itbNzCol + pos;
                itColPtUmf = itbColPtUmf + pos;

                std::transform(itNzCol, iteNzCol, itColPtUmf, itColPtUmf+1, std::plus<Idx>());

                itLRow = listRowIdx.begin() ;
                itLVal = listValues.begin();
                std::advance(itLRow, *itColPtUmf);
                std::advance(itLVal, *itColPtUmf);

                listValues.insert(itLVal, *itval);
                listRowIdx.insert(itLRow, rowIdx);

                ++itval;
            }
        }
    }

    resultUmf.resize(listValues.size());
    std::copy(listValues.begin(), listValues.end(), resultUmf.begin());

    rowIdxUmf.resize(listRowIdx.size());
    std::copy(listRowIdx.begin(), listRowIdx.end(), rowIdxUmf.begin());
}

} // end of namespace xlifepp