largeMatrix/csStorage/CsStorage.hpp

/*
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 CsStorage.hpp
  \authors D. Martin, E. Lunéville, NGUYEN Manh Ha, N. Kielbasiewicz
  \since 21 juin 2011
  \date 21 feb 2013

  \brief Definition of the xlifepp::CsStorage class

  xlifepp::CsStorage class is the abstract mother class of all compressed sparse storage methods for large matrices
  \verbatim
                      | RowCsStorage   (like CSR)
  CsStorage --------> | ColCsStorage   (like CSC)
                      | DualCsStorage  (like CSR/CSC)
                      | SymCsStorage   (like CSR)
  \endverbatim

  Only _sym access type storage can be used for a matrix with a symmetry
  property ( _symmetric, _skewSymmetric, _selfAdjoint or _skewAdjoint)
*/
#ifndef CS_STORAGE_HPP
#define CS_STORAGE_HPP

#include "config.h"
#include "../MatrixStorage.hpp"

#ifdef XLIFEPP_WITH_OMP
#include <omp.h>
#endif

namespace xlifepp
{

/*!
   \class CsStorage
   abstract base class of all compressed sparse storage classes
 */
class CsStorage : public MatrixStorage
{
  public:
    // Constructors, Destructor
    CsStorage(AccessType = _dual, string_t id = "CsStorage");               //!< constructor by access type
    CsStorage(number_t, AccessType = _dual, string_t id = "CsStorage");       //!< constructor by access type, number of columns and rows (same)
    CsStorage(number_t, number_t, AccessType = _dual, string_t id = "CsStorage"); //!< constructor by access type, number of columns and rows
    virtual ~CsStorage() {}                                               //!< virtual destructor

    // public utilities
    // virtual void visual(ostream&, number_t vb = 0) const = 0;    //!< to visalize non zero values
    virtual number_t size() const = 0;                       //!< storage size
    virtual number_t lowerPartSize() const {return 0;}       //!< size of lower triangular part except for Dual/SymCsStorage
    virtual number_t upperPartSize() const {return 0;}       //!< size of upper triangular part except for Dual/SymCsStorage

  protected:
    // build utilities for child class
//    void buildCsStorage(const std::vector<std::vector<number_t> >&,
//                        std::vector<number_t>&, std::vector<number_t>&); //!< build storage vectors of row and col storage
    template<class L>
    void buildCsStorage(const std::vector<L>&, std::vector<number_t>&, std::vector<number_t>&);  //!< build storage vectors of row and col storage
    void addCsSubMatrixIndices(std::vector<number_t>&, std::vector<number_t>&,
                               const std::vector<number_t>&,const std::vector<number_t>&,
                               bool lower=true, bool diag=false);      //!< add in csstorage submatrix indices given by its row and column indices (>= 1)
    void toScalarCs(const std::vector<number_t>&, const std::vector<number_t>&,
                    dimen_t, dimen_t, std::vector<number_t>&, std::vector<number_t>&); //!< go to scalar cs storage from non scalar cs storage

    MatrixStorage* toDual();            //!< create a new dual cs storage from current sym cs storage
    CsStorage* toColStorage();          //!< create a new col cs storage from current cs storage
    CsStorage* toRowStorage();          //!< create a new row cs storage from current cs storage


    //! remove rows from matrix
    template<typename T>
    void deleteRowsT(std::vector<number_t>&, std::vector<number_t>&, number_t&, number_t&, number_t, number_t, std::vector<T>&);
    //! remove cols from matrix
    template<typename T>
    void deleteColsT(std::vector<number_t>&, std::vector<number_t>&, number_t&, number_t&, number_t, number_t, std::vector<T>&);

    //skyline conversion utilities
    template<typename Iterator, typename IteratorS>
    void fillSkylineDiagonalPart(Iterator&, IteratorS&) const;  //!< copy diagonal part
    template<typename Iterator, typename IteratorS>
    void fillSkylineTriangularPart(const std::vector<number_t>&, const std::vector<number_t>&,
                                   Iterator&, IteratorS&) const; //!< copy triangular (lower or upper) part as a skyline storage

    // print utilites for child class

    template<typename Iterator>
    void printEntriesAll(StrucType, Iterator&, const std::vector<number_t>&, const std::vector<number_t>&,
                         number_t, number_t, number_t, const string_t&, number_t, std::ostream&) const; //!< print matrix in row/col compressed storage
    template<typename Iterator>
    void printEntriesTriangularPart(StrucType, Iterator&, Iterator&, const std::vector<number_t>&,
                                    const std::vector<number_t>&, number_t, number_t, number_t,
                                    const string_t&, number_t, std::ostream&) const;              //!< print lower/upper part of matrix in sym/dual compressed storage
    template<typename Iterator>
    void printCooTriangularPart(std::ostream&, Iterator&, const std::vector<number_t>&,
                                const std::vector<number_t>&, bool, SymType sym = _noSymmetry) const; //!< print lower/upper part in coordinate format
    template<typename Iterator>
    void printCooDiagonalPart(std::ostream&, Iterator&, number_t) const;                            //!< print diagonal part in coordinate format

    //! load cs matrix from dense file
    template <typename T>
    void loadCsFromFileDense(std::istream&, std::vector<T>&, std::vector<number_t>&,
                             std::vector<number_t>&, SymType, bool);
    //! load cs matrix from coordinate file
    template <typename T>
    void loadCsFromFileCoo(std::istream&, std::vector<T>&, std::vector<number_t>&,
                           std::vector<number_t>&, SymType, bool);
    //! load cs matrix from dense file
    template <typename T>
    void loadCsFromFileDense(std::istream&, std::vector<T>&, std::vector<number_t>&,
                             std::vector<number_t>&, std::vector<number_t>&, std::vector<number_t>&, SymType, bool);
    //! load cs matrix from coordinate file
    template <typename T>
    void loadCsFromFileCoo(std::istream&, std::vector<T>&, std::vector<number_t>&, std::vector<number_t>&,
                           std::vector<number_t>&, std::vector<number_t>&, SymType, bool);

    // product of any part of Cs matrix and vector common to different Cs storages (child classes)
    //! diag matrix * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void diagonalMatrixVector(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const;
    //! vector * diag matrix
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void diagonalVectorMatrix(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const;
    //! lower part * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void lowerMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                           MatIterator& itm, VecIterator& itvb, ResIterator& itrb, SymType sym = _noSymmetry) const;


    //! vector * lower part
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void lowerVectorMatrix(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                           MatIterator& itm, VecIterator& itvb, ResIterator& itrb, SymType sym = _noSymmetry) const;
    //! upper part * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void upperMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                           MatIterator& itm, VecIterator& itvb, ResIterator& itrb, SymType sym = _noSymmetry) const;
    //! vector * upper part
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void upperVectorMatrix(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                           MatIterator& itm, VecIterator& itvb, ResIterator& itrb, SymType sym = _noSymmetry) const;
    //! row matrix * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void rowMatrixVector(const std::vector<number_t>&, const std::vector<number_t>&, MatIterator&, VecIterator&, ResIterator&)  const;
    //! vector * row matrix
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void rowVectorMatrix(const std::vector<number_t>&, const std::vector<number_t>&, MatIterator&, VecIterator&, ResIterator&)  const;
    //! col matrix * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void columnMatrixVector(const std::vector<number_t>&, const std::vector<number_t>&, MatIterator&, VecIterator&, ResIterator&)  const;
    //! vector * col matrix
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void columnVectorMatrix(const std::vector<number_t>&, const std::vector<number_t>&, MatIterator&, VecIterator&, ResIterator&)  const;

    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void bzSorDiagonalMatrixVector(MatIterator& itd, const VecIterator& itvb, ResIterator& itrb, const real_t w) const;

    template<typename MatIterator, typename VecIterator, typename XIterator>
    void bzSorLowerSolver(const MatIterator&, const MatIterator&, VecIterator&, XIterator&, XIterator&, const std::vector<number_t>&, const std::vector<number_t>&, const real_t) const;

    template<typename MatIterator, typename VecIterator, typename XIterator>
    void bzLowerD1Solver(const MatIterator&, const VecIterator&, const XIterator&, const XIterator&, const std::vector<number_t>&, const std::vector<number_t>&) const;

    template<typename MatIterator, typename VecIterator, typename XIterator>
    void bzSorDiagonalSolver(const MatIterator& itdb, VecIterator& itbb, XIterator& itxb, XIterator& itxe, const real_t w) const;

    template<typename MatRevIterator, typename VecRevIterator, typename XRevIterator>
    void bzSorUpperSolver(const MatRevIterator&, const MatRevIterator&, VecRevIterator&, XRevIterator&, XRevIterator&, const std::vector<number_t>&,
                          const std::vector<number_t>&, const real_t, SymType =_noSymmetry) const;

    template<typename MatRevIterator, typename VecRevIterator, typename XRevIterator>
    void bzSorUpperD1Solver(const MatRevIterator&, const MatRevIterator&, VecRevIterator&, XRevIterator&, XRevIterator&, const std::vector<number_t>&,
                            const std::vector<number_t>&, const real_t, SymType =_noSymmetry) const;

    template<typename MatRevIterator, typename VecRevIterator, typename XRevIterator>
    void bzUpperD1Solver(const MatRevIterator&, const MatRevIterator&, VecRevIterator&, XRevIterator&, XRevIterator&, const std::vector<number_t>&,
                         const std::vector<number_t>&, const real_t, SymType =_noSymmetry) const;
    //! sum of any part of dense matrix and dense matrix
    template<typename Mat1Iterator, typename Mat2Iterator, typename ResIterator>
    void sumMatrixMatrix(Mat1Iterator&, Mat2Iterator&, ResIterator&, ResIterator&) const;

    //! extract index thread for block of column or row
    void extractThreadIndex(const std::vector<number_t>& pointer, const std::vector<number_t>& index,
                            number_t& numThread,
                            std::vector<std::vector<number_t>::const_iterator>& itLowerThread,
                            std::vector<std::vector<number_t>::const_iterator>& itUpperThread) const;

#ifdef XLIFEPP_WITH_OMP
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void parallelRowMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                 MatIterator& itm, VecIterator& itvb, ResIterator& itrb,
                                 number_t numThread,
                                 std::vector<std::vector<number_t>::const_iterator>& itThreadLower,
                                 std::vector<std::vector<number_t>::const_iterator>& itThreadUpper) const;

    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void parallelColumnMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                    MatIterator& itm, VecIterator& itvb, ResIterator& itrb,
                                    number_t numThread,
                                    std::vector<std::vector<number_t>::const_iterator>& itThreadLower,
                                    std::vector<std::vector<number_t>::const_iterator>& itThreadUpper,
                                    number_t numRes) const;

    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void parallelLowerMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                   MatIterator& itm, VecIterator& itvb, ResIterator& itrb,
                                   number_t numThread,
                                   std::vector<std::vector<number_t>::const_iterator>& itThreadLower,
                                   std::vector<std::vector<number_t>::const_iterator>& itThreadUpper,
                                   SymType sym) const;

    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void parallelUpperMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                   MatIterator& itm, VecIterator& itvb, ResIterator& itrb,
                                   number_t numThread,
                                   std::vector<std::vector<number_t>::const_iterator>& itThreadLower,
                                   std::vector<std::vector<number_t>::const_iterator>& itThreadUpper,
                                   number_t numRes,
                                   SymType sym) const;
#endif
};


std::vector<number_t> skylinePointer(const std::vector<number_t>&, const std::vector<number_t>&); //!< get the skyline pointers from cs pointers

inline
bool smallPivot(real_t piv)
{
  return (piv<= theZeroThreshold);
}
inline
bool smallPivot(complex_t piv)
{
  return (std::abs(piv)<= theZeroThreshold);
}

/*
--------------------------------------------------------------------------------
  Template functions
--------------------------------------------------------------------------------
*/

/*--------------------------------------------------------------------------------
 template constructor from indices list
--------------------------------------------------------------------------------*/
// construct row (resp. col) storage vectors (index and pointer)  from the indices of columns (resp rows)
// stored in a vector<L > (rcIndex), called by CsStorage child class
// L if template class supporting ::const_iterator and size() member function either vector<number_t>, set<number_t>, list<number_t>, ...
//
// note : rcIndex starts at 1 !!!
template<class L>
void CsStorage::buildCsStorage(const std::vector<L>& rcIndex, std::vector<number_t>& index, std::vector<number_t>& pointer)
{
  trace_p->push("CsStorage::buildCsStorage");
  //update row/colPointer_
  typename std::vector<L>::const_iterator itvs;
  number_t nnz = 0;
  pointer.resize(rcIndex.size() + 1,0);    //+1 to store the size of stored values (more convenient)
  std::vector<number_t>::iterator itvn = pointer.begin();
  for(itvs = rcIndex.begin(); itvs != rcIndex.end(); ++itvs, ++itvn)
    {
      *itvn = nnz;
      nnz += itvs->size();
    }
  *itvn = nnz; // last of rowPointers_ is the number of non zeros
  // update col/rowIndex_
  index.resize(nnz);
  itvn = index.begin();
  typename L::const_iterator its;
  for(itvs = rcIndex.begin(); itvs != rcIndex.end(); ++itvs)
    for(its = itvs->begin(); its != itvs->end(); ++its, ++itvn) { *itvn = *its -1; }
  // end of construction
  trace_p->pop();
}

/*--------------------------------------------------------------------------------
 template delete Row/Col operations
--------------------------------------------------------------------------------*/

// template delete rows/cols rc1,..., rc2, works for ColCsStorage and RowCsStorage
template<typename T>
void CsStorage::deleteRowsT(std::vector<number_t>& crPointer, std::vector<number_t>& rcIndex,
                            number_t& nbcr, number_t& nbrc, number_t rc1in, number_t rc2in, std::vector<T>& v)
{
  number_t rc1=std::min(std::max(rc1in,number_t(1)),nbrc);  //to prevent oversized indices
  number_t rc2=std::min(std::max(rc2in,number_t(1)),nbrc);
  if(rc2<rc1) return; //nothing to remove

  //recreate all storage pointer (expansive)
  std::vector<number_t>::iterator itr=rcIndex.begin();
  std::vector<number_t>::iterator itp=crPointer.begin();
  std::vector<std::vector<number_t> > rowcols(nbcr);
  std::vector<std::vector<number_t> >::iterator itrc=rowcols.begin();
  typename std::vector<T>::iterator itv=v.begin()+1;
  typename std::vector<T>::iterator itnv=itv;
  number_t s=rc2-rc1;
  for(number_t c=0; c<nbcr; c++, itp++, itrc++)
    {
      number_t l=*(itp+1) -*itp;
      std::vector<number_t>& colrow=*itrc;
      for(number_t k=0; k<l; k++, itr++, itv++)
        {
          number_t rc = *itr + 1;
          if(rc<rc1)
            {
              colrow.push_back(*itr + 1);
              *itnv = *itv; itnv++;
            }
          if(rc>rc2)
            {
              colrow.push_back(*itr - s);
              *itnv = *itv; itnv++;
            }
        }
    }
  std::vector<number_t> newrcIndex, newcrPointer;
  buildCsStorage(rowcols, newrcIndex, newcrPointer); //construct new rowIndex and colPointer
  //update storage
  crPointer=newcrPointer;
  rcIndex=newrcIndex;
  nbrc -= rc2-rc1+1;
  v.resize(crPointer[nbcr]+1);
}

// template delete cols/rows cr1,..., cr2, works for ColCsStorage and RowCsStorage
template<typename T>
void CsStorage::deleteColsT(std::vector<number_t>& crPointer, std::vector<number_t>& rcIndex,
                            number_t& nbcr, number_t& nbrc, number_t cr1in, number_t cr2in, std::vector<T>& v)
{

  number_t cr1=std::min(std::max(cr1in,number_t(1)),nbcr);  //to prevent oversized indices
  number_t cr2=std::min(std::max(cr2in,number_t(1)),nbcr);
  if(cr2<cr1) return;

  if(cr2==nbcr) //remove last columns/rows (simple case)
    {
      if(cr1<=1)  //remove all ?
        {
          warning("free_warning","CsStorage::deleteColsT : removing all columns of a matrix !?");
          v.resize(1);
          crPointer.resize(1,0);
          rcIndex.clear();
          nbcr=0;
          return;
        }
      //remove last columns from c1 to nbcols
      crPointer.resize(cr1);
      number_t last=crPointer[cr1-1];
      rcIndex.resize(last);
      v.resize(last+1);
      nbcr=cr1-1;
      return;
    }
  //remove internal columns
  number_t nbk=cr2-cr1+1; //nb of cols/rows to remove
  number_t ad1=crPointer[cr1-1], ad2=crPointer[cr2], ade=crPointer[nbcr];
  number_t nbt=ad2-ad1; //nb of terms to remove
  std::vector<number_t>::iterator itr1=rcIndex.begin()+ad1, itr2=rcIndex.begin()+ad2;
  typename std::vector<T>::iterator itv1=v.begin()+ad1+1, itv2=v.begin()+ad2+1;
  for(number_t k=ad2; k<ade; k++, itr1++, itr2++, itv1++, itv2++)
    {
      *itv1 = *itv2;
      *itr1 = *itr2;
    }
  std::vector<number_t>::iterator itcp=crPointer.begin()+cr2+1;
  for(; itcp!=crPointer.end(); itcp++) *(itcp - nbk) = *itcp -nbt;
  nbcr-=nbk;
  crPointer.resize(nbcr+1);
  number_t last=crPointer[nbcr];
  rcIndex.resize(last);
  v.resize(last+1);
}

/*--------------------------------------------------------------------------------
 template fill a skyline storage
--------------------------------------------------------------------------------*/

template<typename Iterator, typename IteratorS>
void  CsStorage::fillSkylineDiagonalPart(Iterator& it, IteratorS& its) const
{
  for(number_t i = 0; i < diagonalSize(); i++, its++, it++) *its = *it;
}

/*! fill a skyline storage (triangular part) from a cs storage
  it is assumed that :
  - column/row indexes are ordered in cs storage
  - the entries of skyline storage are initialized to 0
  Remark : this function does not required the skyline storage, it is built by an other function

  rcPointer : row or column pointers of cs storage
  crIndex : column or row indexes of cs storage
  it  : iterator on cs storage entries
  its : iterator on skyline storage entries
*/

template<typename Iterator, typename IteratorS>
void  CsStorage::fillSkylineTriangularPart(const std::vector<number_t>& rcPointer, const std::vector<number_t>& crIndex,
    Iterator& it, IteratorS& its) const
{
  std::vector<number_t>::const_iterator itr, itc = crIndex.begin();//, ite = crIndex.end();
  number_t i = 0;
  for(itr = rcPointer.begin(); itr != rcPointer.end() - 1; itr++, i++)
    {
      number_t l = *(itr + 1) - *itr;
      if(l > 0)
        {
          number_t cmin = *itc;
          for(number_t k = 0; k < l; k++, itc++, it++)
            {
              *(its + (*itc - cmin)) = *it;
            }
          its += i - cmin;
        }
    }
}

/*--------------------------------------------------------------------------------
 template print utilities
--------------------------------------------------------------------------------*/
// print matrix of scalars for row or cs strorage
template<typename Iterator>
void CsStorage::printEntriesAll(StrucType st, Iterator& itm, const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                number_t perRow, number_t width, number_t prec, const string_t& rowOrcol, number_t vb, std::ostream& os) const
{
  number_t rcmin = std::min(vb, number_t(pointer.size()) - 1); // number of rows or cols to print
  string_t f = "firste";
  if(rcmin > 1) { f = "firstes"; }
  os << "(" << words(f) << " " << rcmin << " " << words(rowOrcol) << "s.)";
  os.setf(std::ios::scientific);
  string_t colOrrow;
  if(rowOrcol == "row") { colOrrow = "col"; }
  else { colOrrow = "row"; }
  for(number_t rc = 0; rc < rcmin; rc++)
    {
      number_t nnz = pointer[rc + 1] - pointer[rc]; // number of non zeros on row/col rc
      os << eol << "  " << words(rowOrcol) << " " << rc + 1 << " (" << nnz;
      if(nnz == 0) { os << words("no entry") << " )"; }
      if(nnz > 1) { os << " " << words("entries") << ", " << words(colOrrow) << " :"; }
      else        { os << " " << words("entry") << ", " << words(colOrrow) << " :"; }
      if(nnz > 0)
        {
          for(number_t i = pointer[rc]; i < pointer[rc + 1]; i++)  { os << " " << index[i] + 1; }
          os << ")";
          if(st == _scalar) { printRowWise(os, "   ", perRow, width, prec, itm, itm + nnz); }
          else for(Iterator it = itm; it < itm + nnz; it++) { os << *it; }
          itm += nnz;
        }
    }
  os.unsetf(std::ios::scientific);
  os << eol;
}

/*! print matrix of scalars for sym or dual cs strorage
    itd   : iterator on matrix diagonal entries
    itlu  : iterator on matrix lower or upper triangular part entries
    index : numbers of cols (resp. rows) on rows (resp. cols)
    pointer : adresses of the beginning of rows (resp. cols)
*/
template<typename Iterator>
void CsStorage::printEntriesTriangularPart(StrucType st, Iterator& itd, Iterator& itlu, const std::vector<number_t>& index, const std::vector<number_t>& pointer, number_t perRow, number_t width, number_t prec, const string_t& rowOrcol, number_t vb, std::ostream& os) const
{
  number_t rcmin = std::min(vb, number_t(pointer.size()) - 1); // number of rows or cols to print
  string_t f = "firste";
  if(rcmin > 1) { f = "firstes"; }
  os << "(" << words(f) << " " << rcmin << " " << words(rowOrcol) << "s.)";
  os.setf(std::ios::scientific);
  string_t colOrrow;
  if(rowOrcol == "row") { colOrrow = "col"; }
  else { colOrrow = "row"; }
  for(number_t rc = 0; rc < rcmin; rc++)
    {
      number_t nnz = pointer[rc + 1] - pointer[rc]; // number of non zeros on row/col rc
      os << eol << "  " << words(rowOrcol) << " " << rc + 1;
      if(nnz == 0)  // only diagonal coef.
        {
          os << " (1 " << words("entry") << ", " << words(colOrrow) << " : " << rc + 1 << ")";
          if(st == _scalar) { printRowWise(os, "   ", perRow - 1, width, prec, itd, itd + 1); }
          else { os << *itd; }
          itd++;
        }
      if(nnz > 0)   // non zero coefs. and diagonal coef.
        {
          os << " (" << nnz + 1 << " " << words("entries") << ", " << words(colOrrow) << " : ";
          for(number_t i = pointer[rc]; i < pointer[rc + 1]; i++)  { os << " " << index[i] + 1; }
          os << " " << (rc + 1) << ")";
          if(st == _scalar)
            {
              printRowWise(os, "   ", perRow - 1, width, prec, itlu, itlu + nnz);
              os << std::setw(width) << std::setprecision(prec) << *itd++;
            }
          else
            {
              for(Iterator it = itlu; it < itlu + nnz; it++) { os << *it; }
              os << *itd++;
            }
          itlu += nnz;
        }
    }
  os.unsetf(std::ios::scientific);
  os << eol;
}

/*!
 print "triangular" part of matrix to ostream in coordinate form  i j A_ij :
    1 1 A11
    1 2 A12
    ...
 a null value or a value less than a tolerance in module is not displayed
 The coordinates are not ordered
*/
template<typename Iterator>
void CsStorage::printCooTriangularPart(std::ostream& os, Iterator& itm, const std::vector<number_t>& index, const std::vector<number_t>& pointer, bool byRow, SymType sym) const
{
  std::vector<number_t>::const_iterator itr = pointer.begin();
  std::vector<number_t>::const_iterator itc = index.begin();
  number_t k, l, nbr = pointer.size() - 1;

  for(number_t i = 1; i <= nbr; i++, itr++)      // loop on rows/cols
    {
      number_t nnz = *(itr + 1) - *itr;
      for(number_t j = 0; j < nnz; j++, itc++, itm++) // loop on col/row indices
        {
          if(byRow) {k = i; l = *itc + 1;}
          else {k = *itc + 1; l = i;}
          switch(sym)
            {
              case _skewSymmetric : printCoo(os, -*itm, k, l, 0.); break;
              case _skewAdjoint   : printCoo(os, -conj(*itm), k, l, 0.); break;
              case _selfAdjoint   : printCoo(os, conj(*itm), k, l, 0.); break;
              default             : printCoo(os, *itm, k, l, 0.);
            }
        }
    }
}

template<typename Iterator>
void CsStorage::printCooDiagonalPart(std::ostream& os, Iterator& itm, number_t n) const
{
  for(number_t i = 1; i <= n; i++, itm++)   { printCoo(os, *itm, i, i, 0.); }
}


/*--------------------------------------------------------------------------------
 template load from file and save to file utilities
--------------------------------------------------------------------------------*/

/*! load matrix from file including a dense matrix
    A11 A12  ...                            Re(A11) Im(A11) Re(A12) IM(A12) ....
    A21 A22  ...    or if complex values    Re(A21) Im(A21) Re(A22) IM(A22) ....
    ...                                     ...
 values in module less than the tolerance are considered to be outside the storage
 note : algorithms use pos(i,j) function; they are not optimal!
 arguments :
   ifs       : istream for reading
   mat       : values of matrix
   crIndex   : col or row index of non zero coefficients
   rcPointer : row or col pointer in crIndex
   sym       : symmetry property of the matrix (_noSymmetry, _symmetric, ...)
   realAsCmplx : if true, read a real and cast it to complex
 call only for row or col compressed storage
*/

template <typename T>
void CsStorage::loadCsFromFileDense(std::istream& ifs, std::vector<T>& mat, std::vector<number_t>& crIndex, std::vector<number_t>& rcPointer, SymType sym, bool realAsCmplx)
{
  trace_p->push("CsStorage::loadCsFromFileDense");
  if(accessType_ != _row && accessType_ != _col) { error("storage_not_handled", words("storage type",_cs), words("access type",accessType_)); }

  // read the file and construct the colomn or row indices of "non zeros"
  std::vector<std::vector<number_t> > index;
  std::vector<std::vector<T> > values;
  if(accessType_ == _row)
    {index.resize(nbRows_); values.resize(nbRows_);}
  else
    {index.resize(nbCols_); values.resize(nbCols_);}

  for(number_t i = 0; i < nbRows_; i++)
    for(number_t j = 0; j < nbCols_; j++)
      {
        T v; readItem(ifs, v, realAsCmplx);
        if(std::abs(v) > theTolerance)
          {
            if(accessType_ == _row)
              {index[i].push_back(j+1); values[i].push_back(v);}
            else
              {index[j].push_back(i+1); values[j].push_back(v);}
          }
      }
  // construct compressed storage vectors
  buildCsStorage(index, crIndex, rcPointer);
  mat.resize(size() + 1);
  // assign values in "matrix"
  typename std::vector<T>::iterator itv;
  std::vector<number_t>::iterator iti;
  for(number_t i = 0; i < nbRows_; i++)
    {
      itv = values[i].begin();
      for(iti = index[i].begin(); iti != index[i].end(); iti++, itv++)
        if(sym == _noSymmetry || (sym != _noSymmetry && i >= *iti))
          {
            if(accessType_ == _row) { mat[pos(i + 1, *iti)] = *itv; }
            else { mat[pos(*iti, i+1)] = *itv; }
          }
    }
  trace_p->pop();
}

/*!
 load matrix from file including a matrix in coordinate format
    i j Aij      or if complex values    i j Re(Aij) Im(Aij)
    ...                                     ...
 contrary to dense matrix file, zero values are kept in compressed storage
 note : algorithms use pos(i,j) function; they are not optimal!
 */
template <typename T>
void CsStorage::loadCsFromFileCoo(std::istream& ifs, std::vector<T>& mat, std::vector<number_t>& crIndex, std::vector<number_t>& rcPointer, SymType sym, bool realAsCmplx)
{
  trace_p->push("CsStorage::loadCsFromFileCoo");
  if(accessType_ != _row && accessType_ != _col) { error("storage_not_handled", words("storage type",_cs), words("access type",accessType_)); }

  number_t i, j, nbRows_ = 0, nbCols_ = 0;
  std::map<std::pair<number_t, number_t>, T> values;         // to store matrix values with coordinates access
  while(!ifs.eof())
    {
      ifs >> i >> j;
      T v; readItem(ifs, v, realAsCmplx);
      if(i > nbRows_) { nbRows_ = i; }
      if(j > nbCols_) { nbCols_ = j; }
      values[std::pair<number_t, number_t>(i, j)] = v;
    }
  // build storage
  std::vector < std::vector<number_t> > index; // colum/row indices by row/col
  typename std::map<std::pair<number_t, number_t>, T>::iterator itm;
  if(accessType_ == _row)
    {
      index.resize(nbRows_);
      for(itm = values.begin(); itm != values.end(); itm++)
        { index[itm->first.first - 1].push_back(itm->first.second); }
    }
  else
    {
      index.resize(nbCols_);
      for(itm = values.begin(); itm != values.end(); itm++)
        { index[itm->first.second - 1].push_back(itm->first.first); }
    }
  buildCsStorage(index, crIndex, rcPointer);
  mat.resize(size() + 1);
  // assign matrix values
  for(itm = values.begin(); itm != values.end(); itm++)
    {
      i = itm->first.first;
      j = itm->first.second;
      if(sym == _noSymmetry || (sym != _noSymmetry && i >= j)) { mat[pos(i, j)] = itm->second; }
    }
  trace_p->pop();
}

// call only for dual or sym compressed storage
template <typename T>
void CsStorage::loadCsFromFileDense(std::istream& ifs, std::vector<T>& mat, std::vector<number_t>& colIndex, std::vector<number_t>& rowPointer, std::vector<number_t>& rowIndex, std::vector<number_t>& colPointer, SymType sym, bool realAsCmplx)
{
  trace_p->push("CsStorage::loadCsFromFileDense");
  if(accessType_ != _dual && accessType_ != _sym) { error("storage_not_handled", words("storage type",_cs), words("access type",accessType_)); }

  // read the file and construct the column and row indices of "non zeros"
  std::vector<std::vector<number_t> > cIndex, rIndex;
  std::vector<std::vector<T> > lowValues, upValues;
  std::vector<T> diagValues;
  cIndex.resize(nbRows_); lowValues.resize(nbRows_);
  diagValues.resize(diagonalSize());
  bool upper = (accessType_ == _dual) || (accessType_ == _sym && sym == _noSymmetry);
  if(upper) {rIndex.resize(nbCols_); upValues.resize(nbCols_);}

  for(number_t i = 0; i < nbRows_; i++)
    for(number_t j = 0; j < nbCols_; j++)
      {
        T v; readItem(ifs, v, realAsCmplx);
        if(std::abs(v) > theTolerance)
          {
            if(i > j) {cIndex[i].push_back(j+1); lowValues[i].push_back(v);}
            if(i == j) { diagValues[i] = v; }
            if(i < j && upper) {rIndex[j].push_back(i+1); upValues[j].push_back(v);}
          }
      }
  // construct compressed storage vectors
  buildCsStorage(cIndex, colIndex, rowPointer);             // lower triangular part
  if(upper) { buildCsStorage(rIndex, rowIndex, colPointer); }     // upper triangular part
  if(upper) { mat.resize(diagonalSize() + lowerPartSize() + upperPartSize() + 1); }
  else      { mat.resize(diagonalSize() + lowerPartSize() + 1); }

  // load values in "matrix"
  mat[0] = T();
  typename std::vector<T>::iterator itm = mat.begin() + 1;
  typename std::vector<std::vector<T> >::iterator itvv;
  typename std::vector<T>::iterator itv;
  std::vector<std::vector<number_t> >::iterator itrc;
  std::vector<number_t>::iterator iti;
  for(itv = diagValues.begin(); itv != diagValues.end(); itv++, itm++) { *itm = *itv; } // diagonal part
  itvv = lowValues.begin();

  number_t i = 1;
  for(itrc = cIndex.begin(); itrc != cIndex.end(); itrc++, itvv++, i++)      // lower triangular part
    {
      itv = itvv->begin();
      for(iti = itrc->begin(); iti != itrc->end(); iti++, itv++)  { mat[pos(i, *iti)] = *itv; }
    }

  if(upper)  // upper triangular part
    {
      itvv = upValues.begin();
      number_t i = 1;
      for(itrc = rIndex.begin(); itrc != rIndex.end(); itrc++, itvv++, i++)
        {
          itv = itvv->begin();
          for(iti = itrc->begin(); iti != itrc->end(); iti++, itv++)  { mat[pos(*iti, i)] = *itv; }
        }
    }

  trace_p->pop();
}

/*!
 load matrix from file including a matrix in coordinate format
    i j Aij      or if complex values    i j Re(Aij) Im(Aij)
    ...                                     ...
 contrary to dense matrix file, zero values are kept in compressed storage
 note : algorithm use pos(i,j) function; it is not optimal!
*/
template <typename T>
void CsStorage::loadCsFromFileCoo(std::istream& ifs, std::vector<T>& mat, std::vector<number_t>& colIndex, std::vector<number_t>& rowPointer, std::vector<number_t>& rowIndex, std::vector<number_t>& colPointer, SymType sym, bool realAsCmplx)
{
  trace_p->push("CsStorage::loadFromFileCoo");
  if(accessType_ != _dual && accessType_ != _sym) { error("storage_not_handled", words("storage type",_cs), words("access type",accessType_)); }

  // read values and store them in a map indexed by coordinates matrix
  number_t i, j, nbRows_ = 0, nbCols_ = 0;
  std::map<std::pair<number_t, number_t>, T> values;         // to store matrix values with coordinates access
  while(!ifs.eof())
    {
      ifs >> i >> j;
      T v; readItem(ifs, v, realAsCmplx);
      if(i > nbRows_) { nbRows_ = i; }
      if(j > nbCols_) { nbCols_ = j; }
      values[std::pair<number_t, number_t>(i, j)] = v;
    }

  // build storage
  bool upper = (accessType_ == _dual) || (accessType_ == _sym && sym == _noSymmetry);
  std::vector < std::vector<number_t> > cIndex, rIndex; //colum and row indices by row/col
  typename std::map<std::pair<number_t, number_t>, T>::iterator itmv;
  cIndex.resize(nbRows_);
  if(upper) { rIndex.resize(nbCols_); }
  for(itmv = values.begin(); itmv != values.end(); itmv++)
    {
      number_t i = itmv->first.first, j = itmv->first.second;
      if(i > j) { cIndex[i - 1].push_back(j); }
      if(i < j && upper) { rIndex[j - 1].push_back(i); }
    }
  buildCsStorage(cIndex, colIndex, rowPointer);
  if(accessType_ == _dual) { buildCsStorage(rIndex, rowIndex, colPointer); }
  if(upper) { mat.resize(diagonalSize() + lowerPartSize() + upperPartSize() + 1); }
  else      { mat.resize(diagonalSize() + lowerPartSize() + 1); }

  // assign matrix values
  for(itmv = values.begin(); itmv != values.end(); itmv++)
    {
      i = itmv->first.first;
      j = itmv->first.second;
      if((i >= j) || (upper && j > i)) { mat[pos(i, j)] = itmv->second; }
    }

  trace_p->pop();
}

/*--------------------------------------------------------------------------------
   template Matrix x Vector partial multiplications used in child classes
 --------------------------------------------------------------------------------*/
/*!template Matrix x Vector partial multiplications used in child classes
   CAUTION ! We use non-commutative * operations here as we may deal
             with overloaded matrix-vector left or right products
   in functions, iterators have the following meaning
     itd  : iterator on matrix diagonal entries (begin)
     itm  : iterator on matrix entries (lower "triangular" part)
     itvb : iterator on given vector entries (begin)
     itve : iterator on given vector entries (end)
     itrb : iterator on result vector entries (begin)
     itre : iterator on result vector entries (end)
     index : col or row indices of non zeros coefs
     pointer : adresses in index of beginning of rows or cols
   Result vector res *** NEED NOT BE INITIALIZED HERE ***
*/

// partial multiplication by diagonal M*v, matrix need not be square
// should be called first in a product matrix*vector (result vector initialisation)
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::diagonalMatrixVector(MatIterator& itd, VecIterator& itvb, ResIterator& itrb, ResIterator& itre) const
{
#ifdef XLIFEPP_WITH_OMP
  MatIterator itm = itd;
  VecIterator itv = itvb;
  ResIterator itr;
  #pragma omp parallel default(none) \
  private(itm, itv, itr) \
  shared(itvb, itd, itrb, itre)
  {
    number_t distance = 0;
    #pragma omp for nowait
    for(itr = itrb; itr < (itrb + diagonalSize()); ++itr)
      {
        distance = itr - itrb;
        itm = itd  + distance;
        itv = itvb + distance;
        *itr = *itm * *itv;
      }

    #pragma omp for nowait
    for(itr = itrb + diagonalSize(); itr < itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  }
#else
  VecIterator itv = itvb;
  ResIterator itr;
  for(itr = itrb; itr != itrb + diagonalSize(); ++itr, ++itd, ++itv) { *itr = *itd * *itv; }
  for(itr = itrb + diagonalSize(); itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
#endif
}

// partial multiplication by diagonal v*M, matrix need not be square
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::diagonalVectorMatrix(MatIterator& itd, VecIterator& itvb, ResIterator& itrb, ResIterator& itre) const
{
#ifdef XLIFEPP_WITH_OMP
  MatIterator itm = itd;
  VecIterator itv = itvb;
  ResIterator itr;
  #pragma omp parallel default(none) \
  private(itm, itv, itr) \
  shared(itvb, itd, itrb, itre)
  {
    number_t distance = 0;
    #pragma omp for nowait
    for(itr = itrb; itr < (itrb + diagonalSize()); ++itr)
      {
        distance = itr - itrb;
        itm = itd  + distance;
        itv = itvb + distance;
        *itr = *itm * *itv;
      }

    #pragma omp for nowait
    for(itr = itrb + diagonalSize(); itr < itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  }
#else
  VecIterator itv = itvb;
  ResIterator itr;
  for(itr = itrb; itr != itrb + diagonalSize(); ++itr, ++itd, ++itv) { *itr = *itv * *itd; }
  for(itr = itrb + diagonalSize(); itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
#endif
}

// partial multiplication by lower triangular part
// also used in multiplication of upper part of transposed matrix by vector
// should be called after Diagonal multiplication
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::lowerMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                  MatIterator& itm, VecIterator& itvb, ResIterator& itrb, SymType sym) const
{
#ifdef XLIFEPP_WITH_OMP
  // Get the number of thread to execute code
  number_t numThread = numberOfThreads();
  const number_t GRANULARITY = 16;
  numThread *= GRANULARITY;
  // Lower and upper bound position for each thread
  std::vector<std::vector<number_t>::const_iterator> itThreadLower(numThread);
  std::vector<std::vector<number_t>::const_iterator> itThreadUpper(numThread);
  extractThreadIndex(pointer, index, numThread, itThreadLower, itThreadUpper);
  parallelLowerMatrixVector(index, pointer, itm, itvb, itrb, numThread, itThreadLower, itThreadUpper, sym);

#else
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;
  ResIterator itr = itrb;
  switch(sym)
    {
      case _skewSymmetric:
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*itr -= *itm * *(itvb + *iti); ++iti; ++itm;}
          }
        break;
      case _selfAdjoint:
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*itr += conj(*itm) * *(itvb + *iti); ++iti; ++itm;}
          }
        break;
      case _skewAdjoint:
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*itr -= conj(*itm) * *(itvb + *iti); ++iti; ++itm;}
          }
        break;
      default: // symmetric matrix
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*itr += *itm * *(itvb + *iti); ++iti; ++itm;}
          }
    }
#endif
}

template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::lowerVectorMatrix(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                  MatIterator& itm, VecIterator& itvb, ResIterator& itrb, SymType sym) const
{
#ifdef XLIFEPP_WITH_OMP
  // Get the number of thread to execute code
  number_t numThread = numberOfThreads();
  const number_t GRANULARITY = 4;
  numThread *= GRANULARITY;
  // Lower and upper bound position for each thread
  std::vector<std::vector<number_t>::const_iterator> itThreadLower(numThread);
  std::vector<std::vector<number_t>::const_iterator> itThreadUpper(numThread);
  extractThreadIndex(pointer, index, numThread, itThreadLower, itThreadUpper);
  parallelUpperMatrixVector(index, pointer, itm, itvb, itrb, numThread, itThreadLower, itThreadUpper, MatrixStorage::nbCols_, sym);
#else
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;
  VecIterator itv = itvb;
  switch(sym)
    {
      case _skewSymmetric:
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*(itrb + *iti) -= *itv * *itm; ++iti; ++itm;}
          }
        break;
      case _selfAdjoint:
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*(itrb + *iti) += *itv * conj(*itm); ++iti; ++itm;}
          }
        break;
      case _skewAdjoint:
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*(itrb + *iti) -= *itv * conj(*itm); ++iti; ++itm;}
          }
        break;
      default: //symmetric matrix
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*(itrb + *iti) += *itv * *itm; ++iti; ++itm;}
          }
        break;
    }
#endif
}

// partial multiplication by upper triangular part
// also used in multiplication of upper part of tranposed matrix by vector
// should be called after Diagonal multiplication
// note : for block matrix, (Matrix<Matrix>) a symmetry property means a symmetry property for the full matrix not for the block matrix
//        thus the product of the upper part has to be "transposed" because Aji=tAij
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::upperMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                  MatIterator& itm, VecIterator& itvb, ResIterator& itrb, SymType sym) const
{
#ifdef XLIFEPP_WITH_OMP
  // Get the number of thread to execute code
  number_t numThread = numberOfThreads();
  const number_t GRANULARITY = 4 ;//16;
  numThread *= GRANULARITY;
  // Lower and upper bound position for each thread
  std::vector<std::vector<number_t>::const_iterator> itThreadLower(numThread);
  std::vector<std::vector<number_t>::const_iterator> itThreadUpper(numThread);
  extractThreadIndex(pointer, index, numThread, itThreadLower, itThreadUpper);
  parallelUpperMatrixVector(index, pointer, itm, itvb, itrb, numThread, itThreadLower, itThreadUpper, MatrixStorage::nbRows_, sym);
#else
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;
  VecIterator itv = itvb;

  switch(sym)
    {
      case _skewSymmetric:
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            //while(iti != itie) {*(itrb + *iti) -= *itm * *itv; iti++; itm++;}
            while(iti != itie) {*(itrb + *iti) -= *itv * *itm ; ++iti; ++itm;} //modif to deal with block matrix which is skew-symmetric, not block skew-symmetric
          }
        break;
      case _selfAdjoint:
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            //while(iti != itie) {*(itrb + *iti) += conj(*itm) * *itv; iti++; itm++;}
            while(iti != itie) {*(itrb + *iti) += *itv * conj(*itm) ; ++iti; ++itm;} //modif to deal with block matrix which is selfadjoint, not block selfadjoint
          }
        break;
      case _skewAdjoint:
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            //while(iti != itie) {*(itrb + *iti) -= conj(*itm) * *itv; iti++; itm++;}
            while(iti != itie) {*(itrb + *iti) -= *itv * conj(*itm) ; ++iti; ++itm;} //modif to deal with block matrix which is skewadjoint, not block skewadjointwhile

          }
        break;
      case _symmetric :
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            //while(iti != itie) { *(itrb + *iti) += *itm * *itv; iti++; itm++;}
            while(iti != itie) { *(itrb + *iti) +=  *itv * *itm; ++iti; ++itm;}  //modif to deal with block matrix which is symmetric, not block symmetric
          }
        break;
      default: // no symmetry
        for(itp = itpb; itp != itpe; ++itp, ++itv)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) { *(itrb + *iti) += *itm * *itv; ++iti; ++itm;}

          }
    }
#endif
}

// note : for block matrix, (Matrix<Matrix>) a symmetry property means a symmetry property for the full matrix not for the block matrix
//        thus the product of the upper part has to be "transposed" because Aji=tAij
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::upperVectorMatrix(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                  MatIterator& itm, VecIterator& itvb, ResIterator& itrb, SymType sym) const
{
#ifdef XLIFEPP_WITH_OMP
  // Get the number of thread to execute code
  number_t numThread = numberOfThreads();
  number_t GRANULARITY = 4;
  numThread *= GRANULARITY;
  // Lower and upper bound position for each thread
  std::vector<std::vector<number_t>::const_iterator> itThreadLower(numThread);
  std::vector<std::vector<number_t>::const_iterator> itThreadUpper(numThread);
  extractThreadIndex(pointer, index, numThread, itThreadLower, itThreadUpper);
  parallelLowerMatrixVector(index, pointer, itm, itvb, itrb, numThread, itThreadLower, itThreadUpper, sym);

#else
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;
  ResIterator itr = itrb;
  switch(sym)
    {
      case _skewSymmetric:
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            //while(iti != itie) {*itr -= *(itvb + *iti) * *itm ; iti++; itm++;}
            while(iti != itie) {*itr -= *itm * *(itvb + *iti) ; ++iti; ++itm;}
          }
        break;
      case _selfAdjoint:
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            //while(iti != itie) {*itr += *(itvb + *iti) * conj(*itm); iti++; itm++;}
            while(iti != itie) {*itr += conj(*itm) * *(itvb + *iti); ++iti; ++itm;}
          }
        break;
      case _skewAdjoint:
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            //while(iti != itie) {*itr -= *(itvb + *iti) * conj(*itm); iti++; itm++;}
            while(iti != itie) {*itr -= conj(*itm) * *(itvb + *iti); ++iti; ++itm;}
          }
        break;
      case _symmetric:
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            //while(iti != itie) {*itr += *(itvb + *iti) * *itm; iti++; itm++;}
            while(iti != itie) {*itr += *itm * *(itvb + *iti); ++iti; ++itm;}
          }
        break;
      default: // no symmetry
        for(itp = itpb; itp != itpe; ++itp, ++itr)
          {
            itie = index.begin() + *(itp + 1);
            while(iti != itie) {*itr += *(itvb + *iti) * *itm; ++iti; ++itm;}
          }
    }
#endif
}

// multiplication M x V of Cs row-major access matrix M and vector V
// also used in multiplication of column-major transposed matrix by vector
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::rowMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer, MatIterator& itm, VecIterator& itvb, ResIterator& itrb)  const
{
#ifdef XLIFEPP_WITH_OMP
  // Get the number of thread to execute code
  number_t numThread = numberOfThreads();
  const number_t GRANULARITY = 16;
  numThread *= GRANULARITY;
  // Lower and upper bound position for each thread
  std::vector<std::vector<number_t>::const_iterator> itThreadLower(numThread);
  std::vector<std::vector<number_t>::const_iterator> itThreadUpper(numThread);
  extractThreadIndex(pointer, index, numThread, itThreadLower, itThreadUpper);
  parallelRowMatrixVector(index, pointer, itm, itvb, itrb, numThread, itThreadLower, itThreadUpper);

#else
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;
  ResIterator itr = itrb;
  for(itp = itpb; itp != itpe; ++itr, ++itp)
    {
      *itr *= 0;                      // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
      itie = index.begin() + *(itp + 1);
      while(iti != itie) {*itr += *itm * *(itvb + *iti); ++iti; ++itm;}
    }
#endif
}

// multiplication V x M of vector V and Cs row-major access matrix M
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::rowVectorMatrix(const std::vector<number_t>& index, const std::vector<number_t>& pointer, MatIterator& itm, VecIterator& itvb, ResIterator& itrb)  const
{
#ifdef XLIFEPP_WITH_OMP
  // Get the number of thread to execute code
  number_t numThread = numberOfThreads();
  const number_t GRANULARITY = 4;
  numThread *= GRANULARITY;

  // Lower and upper bound position for each thread
  std::vector<std::vector<number_t>::const_iterator> itThreadLower(numThread);
  std::vector<std::vector<number_t>::const_iterator> itThreadUpper(numThread);
  extractThreadIndex(pointer, index, numThread, itThreadLower, itThreadUpper);
  parallelColumnMatrixVector(index, pointer, itm, itvb, itrb, numThread, itThreadLower, itThreadUpper, MatrixStorage::nbCols_);

#else
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;
  ResIterator itr = itrb;
  for(itp = itpb; itp != itpe; ++itr, ++itp) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  number_t r = 0;
  itr = itrb;
  for(itp = itpb; itp != itpe; ++itr, ++r, ++itp)
    {
      itie = index.begin() + *(itp + 1);
      while(iti != itie) {*(itrb + *iti) += *(itvb + r) * *itm ; ++iti; ++itm;}
    }
#endif
}

// multiplication M x V of Cs column-major access matrix M and vector V
// also used in multiplication of row-major transposed matrix by vector
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::columnMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer, MatIterator& itm, VecIterator& itvb, ResIterator& itrb) const
{
#ifdef XLIFEPP_WITH_OMP
  // Get the number of thread to execute code
  number_t numThread = numberOfThreads();
  const number_t GRANULARITY = 4 ;//16;
  numThread *= GRANULARITY;
  // Lower and upper bound position for each thread
  std::vector<std::vector<number_t>::const_iterator> itThreadLower(numThread);
  std::vector<std::vector<number_t>::const_iterator> itThreadUpper(numThread);
  extractThreadIndex(pointer, index, numThread, itThreadLower, itThreadUpper);
  parallelColumnMatrixVector(index, pointer, itm, itvb, itrb, numThread, itThreadLower, itThreadUpper, MatrixStorage::nbRows_);

#else

  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;
  VecIterator itv = itvb;
  ResIterator itr = itrb;
  for(itp = itpb; itp != itpe; ++itr, ++itp) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  itr = itrb;
  for(itp = itpb; itp != itpe; ++itv, ++itp)
    {
      itie = index.begin() + *(itp + 1);
      while(iti != itie) {*(itrb + *iti) += *itm * *itv; ++iti; ++itm;}
    }
#endif
}


// multiplication V x M of vector V and Cs column-major access matrix M
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::columnVectorMatrix(const std::vector<number_t>& index, const std::vector<number_t>& pointer, MatIterator& itm, VecIterator& itvb, ResIterator& itrb) const
{
#ifdef XLIFEPP_WITH_OMP
  // Get the number of thread to execute code
  number_t numThread = numberOfThreads();
  const number_t GRANULARITY = 16;
  numThread *= GRANULARITY;
  // Lower and upper bound position for each thread
  std::vector<std::vector<number_t>::const_iterator> itThreadLower(numThread);
  std::vector<std::vector<number_t>::const_iterator> itThreadUpper(numThread);
  extractThreadIndex(pointer, index, numThread, itThreadLower, itThreadUpper);
  parallelRowMatrixVector(index, pointer, itm, itvb, itrb, numThread, itThreadLower, itThreadUpper);

#else
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;
  ResIterator itr = itrb;
  for(itp = itpb; itp != itpe; ++itr, ++itp)
    {
      *itr *= 0;                     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
      itie = index.begin() + *(itp + 1);
      while(iti != itie) {*itr += *(itvb + *iti) * *itm ; ++iti; ++itm;}
    }
#endif
}

/*!
   special template solvers (used for SSOR solvers for Symm & Dual Cs Storages)
    *** Matrix is assumed INVERTIBLE here ***, no check on division by zero is done
    *** SOR-like w*D v diagonal matrix x vector multiplication
*/
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::bzSorDiagonalMatrixVector(MatIterator& itd, const VecIterator& itvb, ResIterator& itrb, const real_t w) const
{
  VecIterator itv = itvb;
  for(ResIterator itr = itrb; itr != itrb + diagonalSize(); ++itr, ++itd, ++itv)
    { *itr = w * (*itd * *itv); }
}

/*!
   SOR lower triangular part (D/w+L) x = b solver
   \param itdb is a (forward) iterator on matrix diagonal entries (D)
   \param itlb is a (forward) iterator on matrix strict lower triangular part entries (L)
   \param itbb is a (forward) iterator on right-hand side vector entries (b)
   \param itxb,itxe is a (forward) begin (resp. end) iterator on solution vector entries (x)
   \param index
   \param pointer
   \param w w coefficient
*/
template<typename MatIterator, typename VecIterator, typename XIterator>
void CsStorage::bzSorLowerSolver(const MatIterator& itdb, const MatIterator& itlb, VecIterator& itbb, XIterator& itxb, XIterator& itxe, const std::vector<number_t>& index, const std::vector<number_t>& pointer, const real_t w) const
{
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;

  MatIterator itd=itdb, itl=itlb;
  VecIterator itb = itbb;
  XIterator itx = itxb;

  for(itp = itpb; itp != itpe; ++itp, ++itx)
    {
      *itx = *itb++;

      itie = index.begin() + *(itp + 1);
      while(iti != itie) {*itx -= *itl * *(itxb + *iti); ++iti; ++itl; }
      *itx *= w / *itd++;
    }
}

/*!
   lower triangular part L x = b solver with diagonal of L=1
   \param itlb is a (forward) iterator on matrix strict lower triangular part entries (L)
   \param itbb is a (forward) iterator on right-hand side vector entries (b)
   \param itxb,itxe is a (forward) begin (resp. end) iterator on solution vector entries (x)
   \param index
   \param pointer
*/
template<typename MatIterator, typename VecIterator, typename XIterator>
void CsStorage::bzLowerD1Solver(const MatIterator& itlb, const VecIterator& itbb, const XIterator& itxb, const XIterator& itxe, const std::vector<number_t>& index, const std::vector<number_t>& pointer) const
{
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end() - 1, itp;
  std::vector<number_t>::const_iterator iti = index.begin(), itie;

  MatIterator itl=itlb;
  VecIterator itb = itbb;
  XIterator itx = itxb;

  for(itp = itpb; itp != itpe; ++itp, ++itx)
    {
      *itx = *itb++;

      itie = index.begin() + *(itp + 1);
      while(iti != itie) {*itx -= *itl * *(itxb + *iti); ++iti; ++itl; }
    }
}

/*!
    SOR-like diagonal D/w x = b solver
    \param itdb is a (forward) iterator on matrix diagonal entries (D)
    \param itbb is a (forward) iterator on right hand side vector b entries
    \param itxb, itxe is a (forward) iterator on solution vector x entries
    \param w w coefficient
*/
template<typename MatIterator, typename VecIterator, typename XIterator>
void CsStorage::bzSorDiagonalSolver(const MatIterator& itdb, VecIterator& itbb, XIterator& itxb, XIterator& itxe, const real_t w) const
{
  MatIterator itd=itdb;
  VecIterator itb=itbb;
  if(w!=1)
    {
      for(XIterator itx = itxb; itx != itxe; ++itx, ++itb, ++itd)
        *itx = w * *itb / *itd;
    }
  else
    {
      for(XIterator itx = itxb; itx != itxe; ++itx, ++itb, ++itd)
        *itx = *itb / *itd;
    }
}

/*!
   SOR upper triangular part (D/w+U) x = b solver
   *** LO AND BEHOLD *** we are using reverse_iterator syntax here ***
   *** loops are thus reversed from end to begin (or rather rbegin to rend)
       even though iterators are "++"ed

   \param itrdb is a (backward) iterator on matrix diagonal entries (D)
     starting on matrix last diagonal entry (m.rbegin)
   \param itrub is a (backward) iterator on matrix strict upper triangular part entries (U)
     starting on matrix last entry
   \param itrbb is a reverse iterator "b.rbegin" iterator on right hand side vector (b) entries
     starting on vector last entry
   \param itrxb,itrxe is a reverse iterator "x.rbegin" (resp. "x.rend") iterator on
     solution vector (x) entries starting on vector last entry
   \param index
   \param pointer
   \param w w coefficient
   \param sym : one of _noSymmetry, _symmetric, _skewSymmetric, _selfAdjoint, _skewAdjoint
*/
template<typename MatRevIterator, typename VecRevIterator, typename XRevIterator>
void CsStorage::bzSorUpperSolver(const MatRevIterator& itrdb, const MatRevIterator& itrub, VecRevIterator& itrbb, XRevIterator& itrxb, XRevIterator& itrxe,
                                 const std::vector<number_t>& index, const std::vector<number_t>& pointer, const real_t w, SymType sym) const
{
  std::vector<number_t>::const_reverse_iterator itrpb = pointer.rbegin(), itrpe = pointer.rend() - 1, itrp;
  std::vector<number_t>::const_reverse_iterator itri = index.rbegin(), itrie = index.rbegin();

  MatRevIterator itrd=itrdb, tiru=itrub;
  VecRevIterator itrb = itrbb;
  XRevIterator itrx = itrxb;

  for(itrx = itrxb; itrx != itrxe; ++itrx) { *itrx = *itrb++; }

  itrx = itrxb;
  for(itrp = itrpb; itrp != itrpe; ++itrp, ++itrx)
    {
      *itrx *= w / *itrd++;
      itrie += *(itrp) - *(itrp + 1);
      switch(sym)
        {
          case _skewSymmetric :
            while(itri != itrie) {*(itrxe - *itri - 1) += *tiru * *(itrx); ++itri; ++tiru; }
            break;
          case _selfAdjoint :
            while(itri != itrie) {*(itrxe - *itri - 1) -= conj(*tiru) * *(itrx); ++itri; ++tiru; }
            break;
          case _skewAdjoint :
            while(itri != itrie) {*(itrxe - *itri - 1) += conj(*tiru) * *(itrx); ++itri; ++tiru; }
            break;
          default :
            while(itri != itrie) {*(itrxe - *itri - 1) -= *tiru * *(itrx); ++itri; ++tiru; }
        }
      itri = itrie;
    }
}

/*!
   SOR upper triangular part I+U) x = b solver
   *** LO AND BEHOLD *** we are using reverse_iterator syntax here ***
   *** loops are thus reversed from end to begin (or rather rbegin to rend)
       even though iterators are "++"ed

   \param itrdb is a (backward) iterator on matrix diagonal entries (D)
     starting on matrix last diagonal entry (m.rbegin)
   \param itrub is a (backward) iterator on matrix strict upper triangular part entries (U)
     starting on matrix last entry
   \param itrbb is a reverse iterator "b.rbegin" iterator on right hand side vector (b) entries
     starting on vector last entry
   \param itrxb,itrxe is a reverse iterator "x.rbegin" (resp. "x.rend") iterator on
     solution vector (x) entries starting on vector last entry
   \param index
   \param pointer
   \param w w coefficient
   \param sym : one of _noSymmetry, _symmetric, _skewSymmetric, _selfAdjoint, _skewAdjoint
*/
template<typename MatRevIterator, typename VecRevIterator, typename XRevIterator>
void CsStorage::bzUpperD1Solver(const MatRevIterator& itrdb, const MatRevIterator& itrub, VecRevIterator& itrbb, XRevIterator& itrxb, XRevIterator& itrxe,  const std::vector<number_t>& index, const std::vector<number_t>& pointer, const real_t w, SymType sym) const
{
  std::vector<number_t>::const_reverse_iterator itrpb = pointer.rbegin(), itrpe = pointer.rend() - 1, itrp;
  std::vector<number_t>::const_reverse_iterator itri = index.rbegin(), itrie = index.rbegin();

  MatRevIterator tiru=itrub;
  VecRevIterator itrb = itrbb;
  XRevIterator itrx = itrxb;

  for(itrx = itrxb; itrx != itrxe; ++itrx) { *itrx = *itrb++; }

  itrx = itrxb;
  for(itrp = itrpb; itrp != itrpe; ++itrp, ++itrx)
    {
      *itrx *=1.;// w / *itrd++;
      itrie += *(itrp) - *(itrp + 1);
      switch(sym)
        {
          case _skewSymmetric :
            while(itri != itrie) {*(itrxe - *itri - 1) += *tiru * *(itrx); ++itri; ++tiru; }
            break;
          case _selfAdjoint :
            while(itri != itrie) {*(itrxe - *itri - 1) -= conj(*tiru) * *(itrx); ++itri; ++tiru; }
            break;
          case _skewAdjoint :
            while(itri != itrie) {*(itrxe - *itri - 1) += conj(*tiru) * *(itrx); ++itri; ++tiru; }
            break;
          default :
            while(itri != itrie) {*(itrxe - *itri - 1) -= *tiru * *(itrx); ++itri; ++tiru; }
        }
      itri = itrie;
    }
}

template<typename Mat1Iterator, typename Mat2Iterator, typename ResIterator>
void CsStorage::sumMatrixMatrix(Mat1Iterator& itm1, Mat2Iterator& itm2, ResIterator& itrb, ResIterator& itre) const
{
  for(ResIterator itr = itrb; itr != itre; ++itr)
    {
      *itr = *itm1 + *itm2;
    }
}

#ifdef XLIFEPP_WITH_OMP
/*!
   Parallelize column matrix-vector multiplication
       (as well as row vector-matrix multiplication)

   \param [in] index vector column index
   \param [in] pointer vector row index
   \param [in] itm iterator pointing to first non zero element in the part
   \param [in] itvb iterator pointing to first element of source vector
   \param [in,out] itrb iterator pointing to first element of vector result
   \param [in] numThread number of thread. It's not necessary to be the same number of OpenMp thread
                   it can be even more depending on the GRANULARITY factor
   \param [in] itThreadLower vector containing lower bound of index for each thread
   \param [in] itThreadUpper vector containing upper bound of index for each thread
*/
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::parallelRowMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
                                        MatIterator& itm, VecIterator& itvb, ResIterator& itrb,
                                        number_t numThread,
                                        std::vector<std::vector<number_t>::const_iterator>& itThreadLower,
                                        std::vector<std::vector<number_t>::const_iterator>& itThreadUpper) const
{
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end()-1;
  std::vector<number_t>::const_iterator itib = index.begin(), itie, iti = itib;
  std::vector<number_t>::const_iterator itp = itpb, itpLower, itpUpper;
  ResIterator itr = itrb;
  MatIterator itim = itm;

  std::vector<std::vector<number_t>::const_iterator>::const_iterator itbLower, itbUpper;
  itbLower = itThreadLower.begin();
  itbUpper = itThreadUpper.begin();
  number_t i = 0;

  #pragma omp parallel for default(none)\
  private(i, itr, iti, itie, itim, itpLower, itpUpper) \
  shared(numThread, itbLower, itbUpper, itpb, itib, itm, itvb, itrb) \
  schedule(dynamic,1)
  for(i = 0; i < numThread; i++)
    {
      itpLower = *(itbLower+i);
      itpUpper = *(itbUpper+i);
      for(; itpLower != itpUpper; ++itpLower)
        {
          itr = itrb + (itpLower - itpb);
          *itr *= 0;
          iti  = itib + *(itpLower);
          itie = itib + *(itpLower + 1);
          itim = itm + *(itpLower);

          while(iti != itie)
            {
              *itr += *(itim) * *(itvb + *iti); ++iti; ++itim;
            }
        }
    }
}

/*!
   Parallelize column matrix-vector multiplication
       (as well as row vector-matrix multiplication)

   \param [in] index vector column index
   \param [in] pointer vector row index
   \param [in] itm iterator pointing to first non zero element in the part
   \param [in] itvb iterator pointing to first element of source vector
   \param [in,out] itrb iterator pointing to first element of vector result
   \param [in] numThread number of thread. It's not necessary to be the same number of OpenMp thread
                   it can be even more depending on the GRANULARITY factor
   \param [in] itThreadLower vector containing lower bound of index for each thread
   \param [in] itThreadUpper vector containing upper bound of index for each thread
   \param [in] nResult size of vector result
*/
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::parallelColumnMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
    MatIterator& itm, VecIterator& itvb, ResIterator& itrb,
    number_t numThread,
    std::vector<std::vector<number_t>::const_iterator>& itThreadLower,
    std::vector<std::vector<number_t>::const_iterator>& itThreadUpper,
    number_t nResult) const
{
  typedef typename IterationVectorTrait<ResIterator>::Type ResType;
  std::vector<number_t>::const_iterator itpb = pointer.begin(), itpe = pointer.end()-1;
  std::vector<number_t>::const_iterator itib = index.begin(), itie, iti = itib;
  std::vector<number_t>::const_iterator itp = itpb, itpLower, itpUpper;
  VecIterator itv = itvb;
  ResIterator itr = itrb;
  MatIterator itim = itm;

  std::vector<std::vector<number_t>::const_iterator>::const_iterator itbLower, itbUpper;
  itbLower = itThreadLower.begin();
  itbUpper = itThreadUpper.begin();

  number_t i = 0;
  #pragma omp parallel  default(none) \
  private(i, itr, iti, itie, itim, itpLower, itpUpper, itv) \
  shared(numThread, itbLower, itbUpper, itpb, itib, itm, itvb, itrb, nResult)
  {
    std::vector<ResType> resTemp(nResult,((*itm) * (*itvb)) * 0.);  //use the first product to correctly initialize resTemp in case of vector of vectors
    typename std::vector<ResType>::iterator itbResTemp = resTemp.begin(), iteResTemp = resTemp.end(), itResTemp = itbResTemp;
    #pragma omp for
    for(i = 0; i < nResult; ++i)
      {
        *(itrb + i) *= 0;
      }

    #pragma omp for nowait schedule(dynamic,1)
    for(i = 0; i < numThread; ++i)
      {
        itpLower = *(itbLower+i);
        itpUpper = *(itbUpper+i);
        for(; itpLower != itpUpper; itpLower++)
          {
            itv = itvb + (itpLower - itpb);
            iti  = itib + *(itpLower);
            itie = itib + *(itpLower + 1);
            itim = itm + *(itpLower);
            while(iti != itie)
              {
                itResTemp   = itbResTemp + (*iti);
                (*itResTemp) += (*itim) * (*itv);
                ++iti; ++itim;
              }
          }
      }
    #pragma omp critical (updateResult)
    {
      for(itResTemp = itbResTemp; itResTemp != iteResTemp; ++itResTemp)
        {
          itr = (itrb + (itResTemp - itbResTemp));
          *itr += *itResTemp;
        }
    }
  }
}

/*!
   Parallelize triangular lower part of matrix-vector multiplication
       (as well as triangular upper part of vector-matrix multiplication)
       The switch-case statement inside the parallel region can decrease the performance
       a little bit; however, in this way, it helps avoiding redundant codes

   \param [in] index vector column index
   \param [in] pointer vector row index
   \param [in] itm iterator pointing to first non zero element in the part
   \param [in] itvb iterator pointing to first element of source vector
   \param [in,out] itrb iterator pointing to first element of vector result
   \param [in] numThread number of thread. It's not necessary to be the same number of OpenMp thread
                   it can be even more depending on the GRANULARITY factor
   \param [in] itThreadLower vector containing lower bound of index for each thread
   \param [in] itThreadUpper vector containing upper bound of index for each thread
   \param [in] sym symmetric type of matrix
*/
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::parallelLowerMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
    MatIterator& itm, VecIterator& itvb, ResIterator& itrb,
    number_t numThread,
    std::vector<std::vector<number_t>::const_iterator>& itThreadLower,
    std::vector<std::vector<number_t>::const_iterator>& itThreadUpper,
    SymType sym) const
{
  std::vector<number_t>::const_iterator itpb = pointer.begin();
  std::vector<number_t>::const_iterator itib = index.begin(), itie, iti = itib;
  std::vector<number_t>::const_iterator itpLower, itpUpper;
  ResIterator itr = itrb;
  MatIterator itim = itm;

  std::vector<std::vector<number_t>::const_iterator>::const_iterator itbLower, itbUpper;
  itbLower = itThreadLower.begin();
  itbUpper = itThreadUpper.begin();

  number_t i = 0;
  #pragma omp parallel default(none)\
  private(i, itr, iti, itie, itim, itpLower, itpUpper) \
  shared(numThread, itbLower, itbUpper, itpb, itib, itm, itvb, itrb, sym)
  {
    switch(sym)
      {
        case _skewSymmetric:
          #pragma omp for schedule(dynamic,1)
          for(i = 0; i < numThread; ++i)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; itpLower++)
                {
                  itr = itrb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      *itr -= *(itim) * *(itvb + *iti); ++iti; ++itim;
                    }
                }
            }
          break;

        case _selfAdjoint:
          #pragma omp for schedule(dynamic,1)
          for(i = 0; i < numThread; i++)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; ++itpLower)
                {
                  itr = itrb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      *itr += conj(*itim) * *(itvb + *iti); ++iti; ++itim;
                    }
                }
            }
          break;

        case _skewAdjoint:
          #pragma omp for schedule(dynamic,1)
          for(i = 0; i < numThread; ++i)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; itpLower++)
                {
                  itr = itrb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      *itr -= conj(*itim) * *(itvb + *iti); ++iti; ++itim;
                    }
                }
            }
          break;

        default: // symmetric matrix
          #pragma omp for schedule(dynamic,1)
          for(i = 0; i < numThread; i++)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; ++itpLower)
                {
                  itr = itrb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      *itr += *(itim) * *(itvb + *iti); ++iti; ++itim;
                    }
                }
            }
          break;
      }
  }
}

/*!
   Parallelize matrix-vector multiplication of triangular upper part
       (as well as vector-matrix multiplication of triangular lower part)
       The switch-case statement inside the parallel region can decrease the performance
       a little bit; however, in this way, it helps avoiding redundant codes

   \param [in] index vector row index
   \param [in] pointer vector column index
   \param [in] itm iterator pointing to first non zero element in the part
   \param [in] itvb iterator pointing to first element of source vector
   \param [in,out] itrb iterator pointing to first element of vector result
   \param [in] numThread number of thread. It's not necessary to be the same number of OpenMp thread
                   it can be even more depending on the GRANULARITY factor
   \param [in] itThreadLower vector containing lower bound of index for each thread
   \param [in] itThreadUpper vector containing upper bound of index for each thread
   \param [in] numRes size of vector result
   \param [in] sym symmetric type of matrix
*/
template<typename MatIterator, typename VecIterator, typename ResIterator>
void CsStorage::parallelUpperMatrixVector(const std::vector<number_t>& index, const std::vector<number_t>& pointer,
    MatIterator& itm, VecIterator& itvb, ResIterator& itrb,
    number_t numThread,
    std::vector<std::vector<number_t>::const_iterator>& itThreadLower,
    std::vector<std::vector<number_t>::const_iterator>& itThreadUpper,
    number_t numRes,
    SymType sym) const
{
  typedef typename IterationVectorTrait<ResIterator>::Type ResType;

  std::vector<number_t>::const_iterator itpb = pointer.begin();
  std::vector<number_t>::const_iterator itib = index.begin(), itie, iti = itib;
  std::vector<number_t>::const_iterator itpLower, itpUpper;
  VecIterator itv = itvb;
  ResIterator itr = itrb;
  MatIterator itim = itm;

  std::vector<std::vector<number_t>::const_iterator>::const_iterator itbLower, itbUpper;
  itbLower = itThreadLower.begin();
  itbUpper = itThreadUpper.begin();

  number_t i = 0;
  #pragma omp parallel  default(none) \
  private(i, itr, iti, itie, itim, itpLower, itpUpper, itv) \
  shared(numThread, itbLower, itbUpper, itpb, itib, itm, itvb, itrb, numRes, sym)
  {
    std::vector<ResType> resTemp(numRes,*itrb * 0.);
    typename std::vector<ResType>::iterator itbResTemp = resTemp.begin(), iteResTemp = resTemp.end(), itResTemp = itbResTemp;

    switch(sym)
      {
        case _skewSymmetric:
          #pragma omp for nowait schedule(dynamic,1)
          for(i = 0; i < numThread; i++)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; ++itpLower)
                {
                  itv = itvb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      itResTemp   = itbResTemp + (*iti);
                      //*itResTemp -= *itim * *(itv);
                      *itResTemp -= *itv * *itim;
                      ++iti; ++itim;
                    }
                }
            }
          break;

        case _selfAdjoint:
          #pragma omp for nowait schedule(dynamic,1)
          for(i = 0; i < numThread; i++)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; ++itpLower)
                {
                  itv = itvb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      itResTemp   = itbResTemp + (*iti);
                      //*itResTemp += conj(*itim) * *(itv);
                      *itResTemp += *itv * conj(*itim);
                      ++iti; ++itim;
                    }
                }
            }
          break;

        case _skewAdjoint:
          #pragma omp for nowait schedule(dynamic,1)
          for(i = 0; i < numThread; i++)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; ++itpLower)
                {
                  itv = itvb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      itResTemp   = itbResTemp + (*iti);
                      //*itResTemp -= conj(*itim) * *(itv);
                      *itResTemp -= *itv * conj(*itim);
                      ++iti; ++itim;
                    }
                }
            }
          break;

        case _symmetric :
          #pragma omp for nowait schedule(dynamic,1)
          for(i = 0; i < numThread; i++)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; ++itpLower)
                {
                  itv = itvb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      itResTemp   = itbResTemp + (*iti);
                      //*itResTemp += *itim * *(itv);
                      *itResTemp += *itv * *itim ;
                      ++iti; ++itim;
                    }
                }
            }
          break;

        default: // no symmetry
          #pragma omp for nowait schedule(dynamic,1)
          for(i = 0; i < numThread; i++)
            {
              itpLower = *(itbLower+i);
              itpUpper = *(itbUpper+i);
              for(; itpLower != itpUpper; ++itpLower)
                {
                  itv = itvb + (itpLower - itpb);
                  iti  = itib + *(itpLower);
                  itie = itib + *(itpLower + 1);
                  itim = itm + *(itpLower);

                  while(iti != itie)
                    {
                      itResTemp   = itbResTemp + (*iti);
                      *itResTemp += *itim * *(itv);
                      ++iti; ++itim;
                    }
                }
            }
          break;

      }
    #pragma omp critical (updateResult)
    {
      for(itResTemp = itbResTemp; itResTemp != iteResTemp; ++itResTemp)
        {
          itr = (itrb + (itResTemp - itbResTemp));
          *itr += *itResTemp;
        }
    }
  }
}
#endif// end of #ifdef XLIFEPP_WITH_OMP

} // end of namespace xlifepp

#endif // CS_STORAGE_HPP