largeMatrix/csStorage/RowCsStorage.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 RowCsStorage.cpp
  \authors D. Martin, E. Lunéville, NGUYEN Manh Ha, N. Kielbasiewicz
  \since 21 jun 2011
  \date 21 feb 2013

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

#include "RowCsStorage.hpp"
#include <functional>

namespace xlifepp
{

/*--------------------------------------------------------------------------------
  Constructors, Destructor
--------------------------------------------------------------------------------*/
//! default constructor
RowCsStorage::RowCsStorage(number_t nr, number_t nc, string_t id)
  : CsStorage(nr, nc, _row, id)
{
  trace_p->push("RowCsStorage constructor");
  std::vector<std::vector<number_t> > cols(nr);
  buildCsStorage(cols, colIndex_, rowPointer_); //construct colIndex_ and rowPointer_
  trace_p->pop();
}

/*! constructor by a pair of numbering vectors rowNumbers and colNumbers (same size)
   rowNumbers[k] is the row numbers of submatrix k
   colNumbers[k] is the col numbers of submatrix k
 indices start from 1 !!!!
 this constructor is well adapted to FE matrix where submatrix is linked to an element
 */
RowCsStorage::RowCsStorage(number_t nr, number_t nc, const std::vector<std::vector<number_t> >& rowNumbers, const std::vector<std::vector<number_t> >& colNumbers, string_t id)
  : CsStorage(nr, nc, _row, id)
{
  trace_p->push("RowCsStorage constructor");
  // build list of col index for each rows
  std::vector<std::set<number_t> > colset(nbRows_);   //temporary vector of col index
  std::vector<std::vector<number_t> >::const_iterator itrs, itcs = colNumbers.begin();
  std::vector<number_t>::const_iterator itr, itc;
  for(itrs = rowNumbers.begin(); itrs != rowNumbers.end(); itrs++, itcs++)
    {
      for(itr = itrs->begin(); itr != itrs->end(); itr++)
        for(itc = itcs->begin(); itc != itcs->end(); itc++)
          { colset[*itr - 1].insert(*itc); }
    }
  // convert set to vector
  std::vector<std::vector<number_t> > cols(nbRows_);   //temporary vector of col index
  std::vector<std::set<number_t> >::iterator itset;
  std::vector<std::vector<number_t> >::iterator itcols = cols.begin();
  for(itset = colset.begin(); itset != colset.end(); itset++, itcols++)
    { *itcols = std::vector<number_t>(itset->begin(), itset->end()); }
  // construct colIndex_ and rowPointer_
  buildCsStorage(cols, colIndex_, rowPointer_);
  trace_p->pop();
}

//! create a new scalar RowCs storage from current RowCs storage and submatrix sizes
RowCsStorage* RowCsStorage::toScalar(dimen_t nbr, dimen_t nbc)
{
  RowCsStorage* cs = new RowCsStorage(nbr*nbRows_, nbc*nbCols_, stringId+"_scalar");
  cs->colIndex_.resize(nbr*nbc*colIndex_.size());
  toScalarCs(rowPointer_, colIndex_, nbc, nbr, cs->rowPointer_, cs->colIndex_);
  return cs;
}

//! check if two storages have the same structures
bool RowCsStorage::sameStorage(const MatrixStorage& sto) const
{
  if(!(sto.storageType()==storageType_ &&
       sto.accessType()==accessType_ &&
       sto.nbOfRows()==nbRows_ &&
       sto.nbOfColumns()==nbCols_&&
       sto.size()==size()))
    return false;
  const RowCsStorage& csto=static_cast<const RowCsStorage&>(sto);
  if(rowPointer_!=csto.rowPointer_) return false;
  if(colIndex_!=csto.colIndex_) return false;
  return true;
}

/*!--------------------------------------------------------------------------------
 add submatrix indices given by its row and column indices (>= 1)
--------------------------------------------------------------------------------*/
void RowCsStorage::addSubMatrixIndices(const std::vector<number_t>& rows,const std::vector<number_t>& cols)
{
  addCsSubMatrixIndices(rowPointer_,colIndex_,rows,cols, false,true);
}

/*!--------------------------------------------------------------------------------
 add in storage indexes (i,j) given by rows (i) and columns (j) (i,j >= 1)
 rows and cols are of same size
--------------------------------------------------------------------------------*/
void RowCsStorage::addIndices(const std::vector<number_t>& rows, const std::vector<number_t>& cols)
{
  if(rows.size()==0) return;  //nothing to do

  std::vector<number_t>::const_iterator itr = rows.begin(),itc=cols.end();
  std::map<number_t, std::set<number_t> > rowsmap;
  std::map<number_t, std::set<number_t> >::iterator itm;;
  //collect col indexes by row
  for(; itr!=rows.end(); itr++, itc++)
    {
      itm=rowsmap.find(*itr);
      if(itm==rowsmap.end()) rowsmap[*itr]=std::set<number_t>();
      rowsmap[*itr].insert(*itc);
    }
  //modify current rowpointer and colindex (expansive operation)
  for(itm=rowsmap.begin(); itm!=rowsmap.end(); itm++)
    if(itm->second.size()>0) addRow(itm->first,itm->second,_all);  //insert col at row r in storage
}

/*!--------------------------------------------------------------------------------
insert in storage indexes (r,c1), (r,c2), .... c1,c2, ... given by cols (r,c>=1)
  expansive operation !
  MatrixPart not managed !
--------------------------------------------------------------------------------*/
void RowCsStorage::addRow(number_t r, const std::set<number_t>& cols, MatrixPart)
{
  if(cols.size()==0) return;
  number_t rb=rowPointer_[r-1], re=rowPointer_[r];
  std::set<number_t> allcols(colIndex_.begin()+rb, colIndex_.begin()+re);
  std::set<number_t>::iterator its=cols.begin();
  for(; its!=cols.end(); its++) allcols.insert(*its -1); //add current cols (-1 shift)
  number_t nc=allcols.size()- (re - rb);
  if(nc==0) return;                      //nothing to do , no new cols

  std::vector<number_t>::iterator itc=colIndex_.begin()+re;
  colIndex_.insert(itc, nc, 0);          //insert nc element initialized to 0
  its=allcols.begin();
  itc=colIndex_.begin()+rb;
  for(; its!=allcols.end(); its++,itc++) *itc=*its; //copy new columns set
  std::vector<number_t>::iterator itr=rowPointer_.begin()+r;
  for(; itr!=rowPointer_.end(); itr++) *itr+=nc; //shift rowPointer
}

/*!--------------------------------------------------------------------------------
 access operator to position of (i,j)>0 coefficients (return 0 if outside)
 This DIRECT access to (i,j) pointer SHOULD NOT be used in matrix computation
--------------------------------------------------------------------------------*/
number_t RowCsStorage::pos(number_t i, number_t j, SymType s) const
{
  if(i == 0 || i > nbRows_ || j == 0 || j > nbCols_) { return 0; }
  for(number_t k = rowPointer_[i - 1]; k < rowPointer_[i]; k++) // find in list of col index of row i
    {
      if(colIndex_[k] == j - 1) { return k + 1; }               // add 1 because values begin at address 1
    }
  return 0;
}

/*!--------------------------------------------------------------------------------
 access to submatrix adresses, useful for assembling matrices
 rows and cols are vector of index (>0) of row and column, no index ordering is assumed !!
 return the vector of adresses (relative number >0 in storage) in a row storage (as Matrix)
 if a row or a col index is out of storage
     an error is raised if errorOn is true (default)
     position is set to 0 else
 if sy!=_noSymmetry only positions of lower part are requested (not used here)
 remark : since the vector storage indexCol_ is ordered for each line, if the indices of columns
          of the submatrix were also, one could write a faster algorithm
--------------------------------------------------------------------------------*/
void RowCsStorage::positions(const std::vector<number_t>& rows, const std::vector<number_t>& cols, std::vector<number_t>& adrs, bool errorOn, SymType sy) const
{
  number_t nba = rows.size() * cols.size();
  if(adrs.size() != nba) { adrs.resize(nba); }
  std::vector<number_t>::iterator itp = adrs.begin();
  std::vector<number_t>::const_iterator itr, itc, itb, ite;
  for(itr = rows.begin(); itr != rows.end(); itr++)  // loop on row index
    {
      number_t i = *itr, rb = rowPointer_[i - 1], re = rowPointer_[i];
      for(itc = cols.begin(); itc != cols.end(); itc++, itp++) // loop on col index
        {
          itb = colIndex_.begin() + rb;
          ite = colIndex_.begin() + re;
          *itp = 0;
          number_t k=0;
          while(itb != ite && *itp == 0) // find in list of col index of row i
            {

              if(*itb == (*itc - 1)) { *itp = rb + k + 1; }
              itb++; k++;
            }
          if(*itp == 0 && errorOn) { error("storage_outofstorage","RowCs", *itr, *itc); }
        }
    }
}

//! get (col indices, adress) of row r in set [c1,c2] ## SymType not used ##
std::vector<std::pair<number_t, number_t> > RowCsStorage::getRow(SymType s,number_t r, number_t c1, number_t c2) const
{
  number_t nbc=c2;
  if(nbc==0) nbc=nbCols_;
  number_t rb=rowPointer_[r-1], re=rowPointer_[r];
  number_t l=re-rb;
  std::vector<std::pair<number_t, number_t> > cols(l);
  std::vector<std::pair<number_t, number_t> >::iterator itcol=cols.begin();
  l=0;
  for(number_t k = rb; k < re; k++) // find in list of col index of row i
    {
      number_t c=colIndex_[k]+1;
      if(c >= c1 && c<=nbc)
        {
          *itcol = std::make_pair(c,k+1);            // add 1 because values begin at address 1
          l++; itcol++;
        }
    }
  cols.resize(l);
  return cols;
}

//! get (row indices, adress) of col c in set [r1,r2] ## SymType not used ##
std::vector<std::pair<number_t, number_t> > RowCsStorage::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();
  number_t l=0;
  for(number_t i=r1; i<=nbr; i++)
    {
      number_t a=pos(i,c);
      if(a!=0)
        {
          *itrow=std::make_pair(i,a);
          l++; itrow++;
        }
    }
  rows.resize(l);
  return rows;
}

//! get col indices of row r in set [c1,c2], specific CsRow algorithm (no getRows, use generic algorithm)
std::set<number_t> RowCsStorage::getCols(number_t r, number_t c1, number_t c2) const
{
  std::set<number_t> colset;
  number_t nbc=c2;
  if(nbc==0) nbc=nbCols_;
  if(nbc<c1) return colset;
  number_t rb=rowPointer_[r-1], re=rowPointer_[r];
  for(number_t k = rb; k < re; k++) // find in list of col index of row i
    {
      number_t c=colIndex_[k]+1;
      if(c >= c1 && c<=nbc) colset.insert(c);
    }
  return colset;
}

//! get col indices of row r in set [c1,c2], specific CsRow algorithm (no getRows, use generic algorithm)
// no allocation, no check (faster)
void RowCsStorage::getColsV(std::vector<number_t>& colsv, number_t& nbcols, number_t r, number_t c1, number_t c2) const
{
  nbcols=0;
  number_t nbc=c2;
  if(nbc==0) nbc=nbCols_;
  if(nbc<c1) return;
  std::vector<number_t>::iterator itv=colsv.begin();
  number_t rb=rowPointer_[r-1], re=rowPointer_[r];
  for(number_t k = rb; k < re; k++) // find in list of col index of row i
    {
      number_t c=colIndex_[k]+1;
      if(c >= c1 && c<=nbc) {*itv=c; itv++; nbcols++;}
    }
}


//! print storage pointers
void RowCsStorage::print(std::ostream& os) const
{
  printHeader(os);
  os<<"row pointer = "<<rowPointer_<<eol;
  os<<"column index = "<<colIndex_<<eol;
}

/*!--------------------------------------------------------------------------------
 print RowCs Storage matrix to output stream (sym is not used here)
--------------------------------------------------------------------------------*/
// print real scalar matrix
void RowCsStorage::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;
  printEntriesAll(_scalar, itm, colIndex_, rowPointer_, entriesPerRow, entryWidth, entryPrec, "row", vb, os);
}

//! print complex scalar matrix
void RowCsStorage::printEntries(std::ostream& os, const std::vector<complex_t>& m, number_t vb, SymType sym) const
{
  std::vector<complex_t>::const_iterator itm = m.begin() + 1;
  printEntriesAll(_scalar, itm, colIndex_, rowPointer_, entriesPerRow / 2, 2 * entryWidth + 1, entryPrec, "row", vb, os);
}

//! large matrix of vectors are not available for the moment, use Matrix n x 1 or 1 x n
void RowCsStorage::printEntries(std::ostream& os, const std::vector< Vector<real_t> >& m, number_t vb, SymType sym) const
{error("storage_novector", "RowCsStorage::printEntries");}

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

//! print matrix of real matrices
void RowCsStorage::printEntries(std::ostream& os, const std::vector< Matrix<real_t> >& m, number_t vb, SymType sym) const
{
  std::vector< Matrix<real_t> >::const_iterator itm = m.begin() + 1;
  printEntriesAll(_matrix, itm, colIndex_, rowPointer_, entriesPerRow, entryWidth, entryPrec, "row", vb, os);
}

//! print matrix of complex matrices
void RowCsStorage::printEntries(std::ostream& os, const std::vector< Matrix<complex_t> >& m, number_t vb, SymType sym) const
{
  std::vector< Matrix<complex_t> >::const_iterator itm = m.begin() + 1;
  printEntriesAll(_matrix, itm, colIndex_, rowPointer_, entriesPerRow / 2, 2 * entryWidth + 1, entryPrec, "row", vb, os);
}

/*!
--------------------------------------------------------------------------------
 print matrix to ostream in coordinate form : i j M_ij
 according to type of values of matrix (overload the generic method in MatrixStorage)
 sym is not used
--------------------------------------------------------------------------------
*/
void RowCsStorage::printCooMatrix(std::ostream& os, const std::vector<real_t>& m, SymType sym) const
{
  std::vector<real_t>::const_iterator itm = m.begin() + 1;
  printCooTriangularPart(os, itm, colIndex_, rowPointer_, true);
}

void RowCsStorage::printCooMatrix(std::ostream& os, const std::vector<complex_t>& m, SymType sym) const
{
  std::vector<complex_t>::const_iterator itm = m.begin() + 1;
  printCooTriangularPart(os, itm, colIndex_, rowPointer_, true);
}

void RowCsStorage::printCooMatrix(std::ostream& os, const std::vector< Matrix<real_t> >& m, SymType sym) const
{
  std::vector< Matrix<real_t> >::const_iterator itm = m.begin() + 1;
  printCooTriangularPart(os, itm, colIndex_, rowPointer_, true);
}

void RowCsStorage::printCooMatrix(std::ostream& os, const std::vector< Matrix<complex_t> >& m, SymType sym) const
{
  std::vector< Matrix<complex_t> >::const_iterator itm = m.begin() + 1;
  printCooTriangularPart(os, itm, colIndex_, rowPointer_, true);
}

/*!--------------------------------------------------------------------------------
 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 RowCsStorage::multMatrixVector(const std::vector<M>& m, V* vp, R* rp) const
{
  trace_p->push("RowCsStorage::multMatrixVector (pointer form)");
  typename std::vector<M>::const_iterator itm = m.begin() + 1;
  rowMatrixVector(colIndex_, rowPointer_, itm, vp, rp);
  trace_p->pop();
}

template<typename M, typename V, typename R>
void RowCsStorage::multMatrixVector(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
  trace_p->push("RowCsStorage::multMatrixVector");
  typename std::vector<M>::const_iterator itm = m.begin() + 1;
  typename std::vector<V>::const_iterator itvb = v.begin();
  typename std::vector<R>::iterator itrb = r.begin();
  rowMatrixVector(colIndex_, rowPointer_, itm, itvb, itrb);
  trace_p->pop();
}

template<typename M, typename V, typename R>
void RowCsStorage::multVectorMatrix(const std::vector<M>& m, V* vp, R* rp) const
{
  trace_p->push("RowCsStorage::multVectorMatrix (pointer form)");
  typename std::vector<M>::const_iterator itm = m.begin() + 1;
  rowMatrixVector(colIndex_, rowPointer_, itm, vp, rp);
  trace_p->pop();
}

template<typename M, typename V, typename R>
void RowCsStorage::multVectorMatrix(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r) const
{
  trace_p->push("RowCsStorage::multVectorMatrix");
  typename std::vector<M>::const_iterator itm = m.begin() + 1;
  typename std::vector<V>::const_iterator itvb = v.begin();
  typename std::vector<R>::iterator itrb = r.begin();
  rowVectorMatrix(colIndex_, rowPointer_, itm, itvb, itrb);
  trace_p->pop();
}

/*!
   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 RowCsStorage::addMatrixMatrix(const std::vector<M1>& m1, const std::vector<M2>& m2, std::vector<R>& r) const
{
  trace_p->push("RowCsStorage::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 itrb = r.begin() + 1, itre = r.end();
  sumMatrixMatrix(itm1, itm2, itrb, itre);
  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 RowCsStorage::toUmfPack(const std::vector<M1>& m1, std::vector<Idx>& colPtUmf, std::vector<Idx>& rowIdxUmf, std::vector<M1>& resultUmf) const
{
  // Reserve memory for resultUmf and set its size to m1.size()-1 (we don't count for 0 in the first position)s
  resultUmf.reserve(m1.size()-1);
  resultUmf.clear();

  // Reserve memory for rowIdxUmf and set its size to 0
  rowIdxUmf.reserve(m1.size()-1);
  rowIdxUmf.clear();

  // Resize column pointer to nbCol+1
  colPtUmf.clear();
  colPtUmf.resize(nbCols_+1);

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

  std::vector<number_t>::const_iterator itbColIdx = colIndex_.begin(), itColIdx = itbColIdx;
  std::vector<number_t>::const_iterator iteColIdx = colIndex_.end();
  std::vector<number_t>::const_iterator itRowPt = rowPointer_.begin();

  typename std::vector<M1>::const_iterator itm = m1.begin()+1;
  typename std::vector<M1>::const_iterator itval;

  for(number_t colIdx = 0; colIdx < nbCols_; ++colIdx,++itColPtUmf)
    {
      itColIdx = itbColIdx;
      Idx nzCol = Idx(0); // Number of non-zero in each "working" column
      while(itColIdx != iteColIdx)
        {
          itColIdx = std::find(itColIdx,iteColIdx, colIdx);
          if(itColIdx != iteColIdx)
            {
              itRowPt = std::find_if(rowPointer_.begin(), rowPointer_.end(),std::bind2nd(std::greater_equal<number_t>(), itColIdx - itbColIdx+1));

              // Update values
              itval = (itColIdx-itbColIdx) + itm;
              resultUmf.push_back(*itval);

              // Update colPtUmf with number of non-zero value
              ++nzCol;

              // Update rowIdxUmf
              rowIdxUmf.push_back(itRowPt-rowPointer_.begin()-1);

              // Increase iterator to make a search in new range
              ++itColIdx;
            }
        }
      // Update colPtUmf with number of non-zero value
      *itColPtUmf += *(itColPtUmf-1) + nzCol;


    }
}
/*!
   incomplete factorization of matrix M stored as matrix F = L U where
   - L is a lower triangular matrix with unit diagonal and is stored as
   the strict lower triangular part of F
   - U is a upper triangular matrix with diagonal stored as diagonal part of F
   and strict upper triangular part stored as the strict upper triangular part of F :

  */
template<typename M>
void RowCsStorage::ilu(std::vector<M>& m, std::vector<M>& f) const
{

  trace_p->push("RowCsStorage::ilu");
  typename std::vector<M>::iterator f_ib=f.begin()+1;
  std::vector<std::pair<number_t, number_t> > Cj;
  std::vector<M> Diag(nbRows_);
  number_t k;
  Diag[0]=*f_ib;
  number_t c,l;

  for (number_t j=0; j < nbRows_; j++) // boucle sur les colonnes
  {
       Cj=getCol(_noSymmetry,j+1,1,nbRows_);
       for (number_t p=1; p < Cj.size(); p++)
         {//Cj[p].first-1 ligne du peme coef non nul  col j
              l=Cj[p].first-1;
              c = rowPointer_[l];
              k = 0;
              while ( (Cj[k].first-1 < l) && (colIndex_[c] < j) )
              {
                  if (colIndex_[c] == Cj[k].first-1)
                       {
                          *(f_ib+Cj[p].second-1) -=*(f_ib+Cj[k].second-1) * (*(f_ib+c));
                           c++;k++;
                        }
                  else if (colIndex_[c] < Cj[k].first-1)   c++;
                  else /*if  (colIndex_[c] > Cj[k].first-1)  */   k++;
              }
            if (Cj[p].first-1 == j)
              {
                 Diag[j]=*(f_ib+Cj[p].second-1);
                 if(std::abs(Diag[j]) < theZeroThreshold) {error("small_pivot") ; }
               }
            else if (Cj[p].first-1 > j) *(f_ib+Cj[p].second-1)/=Diag[j];
         }
   }
  trace_p->pop();
}

} // end of namespace xlifepp