DenseStorage.hpp - OpenGrok cross reference for /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/largeMatrix/denseStorage/DenseStorage.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 DenseStorage.hpp
  \authors D. Martin, E. Lunéville, NGUYEN Manh Ha, N. Kielbasiewicz
  \since 21 jun 2011
  \date 21 feb 2013

  \brief Definition of the xlifepp::DenseStorage class

  xlifepp::DenseStorage class is the abstract mother class of all dense storage methods for large matrices
  \verbatim
                         | RowDenseStorage   (store row[1],row[2] ...)
  DenseStorage --------> | ColDenseStorage   (store col[1],col[2] ...)
                         | DualDenseStorage
                         | SymDenseStorage
  \endverbatim

  Only _symm access type storage can be used for a matrix with a symmetry
  property ( _symmetric, _skewSymmetric, _selfAdjoint or _skewAdjoint)
*/

#ifndef DENSE_STORAGE_HPP
#define DENSE_STORAGE_HPP

#include "config.h"
#include "../MatrixStorage.hpp"
#ifdef XLIFEPP_WITH_OMP
#include "omp.h"
#endif // XLIFEPP_WITH_OMP

namespace xlifepp
{

/*!
   \class DenseStorage
   abstract base class of all dense storage methods
 */

class DenseStorage : public MatrixStorage
{
  public:
    // Constructors, Destructor
    DenseStorage(string_t id="DenseStorage");                                 //!< default constructor
    DenseStorage(AccessType, string_t id="DenseStorage");                     //!< constructor by access type
    DenseStorage(AccessType, number_t, string_t id="DenseStorage");           //!< constructor by access type, number of columns and rows (same)
    DenseStorage(AccessType, number_t, number_t, string_t id="DenseStorage"); //!< constructor by access type, number of columns and rows
    virtual ~DenseStorage() {}                //!< virtual destructor

    MatrixStorage* toDual();                  //!< create a new dual dense storage from sym dense storage

    // utilities
    virtual number_t pos_(number_t i, number_t j, SymType s = _noSymmetry) const   //! overloaded pos fast returns adress of entry (i,j)
    {return pos(i,j,s); }
    virtual number_t size() const {return nbRows_ * nbCols_;} //!< storage size except for SymDenseStorage
    virtual number_t lowerPartSize() const {return 0;}        //!< size of lower triangular part except for Dual/SymDenseStorage
    virtual number_t upperPartSize() const {return 0;}        //!< size of upper triangular part except for Dual/SymDenseStorage
    std::set<number_t> getRows(number_t c, number_t r1=1, number_t r2=0) const;  //!< get row indices of col c in set [r1,r2]
    std::set<number_t> getCols(number_t r, number_t c1=1, number_t c2=0) const;  //!< get col indices of row r in set [c1,c2]
    void getColsV(std::vector<number_t>&, number_t&, number_t r, number_t c1=1, number_t c2=0) const;  //!< get col indices of row r in set [c1,c2]
    void getRowsV(std::vector<number_t>&, number_t&, number_t c, number_t r1=1, number_t r2=0) const;  //!< get row indices of col c in set [r1,r2]
    bool sameStorage(const MatrixStorage&) const;             //!< check if two storages have the same structures

  protected:
    //@{
    //! print utilites for child class
    template<typename Iterator>
    void printScalarEntries(Iterator& it, number_t nrows, number_t ncols, number_t perRow, number_t width, number_t prec,
                            const string_t& rowOrcol, number_t vb, std::ostream& os) const;
    template<typename Iterator>
    void printScalarEntriesTriangularPart(Iterator& itd, Iterator& it, number_t nrows, number_t ncols, number_t perRow, number_t width, number_t prec,
                                          const string_t& rowOrcol, number_t vb, std::ostream& os) const;
    template<typename Iterator>
    void printMatrixEntries(Iterator& it, number_t nrows, number_t ncols,
                            const string_t& rowOrcol, number_t vb, std::ostream& os) const;
    template<typename Iterator>
    void printMatrixEntriesTriangularPart(Iterator& itd, Iterator& it, number_t nrows, number_t ncols,
                                          const string_t& rowOrcol, number_t vb, std::ostream& os) const;
    //@}

    void loadFromFileDense(std::istream&, std::vector<real_t>&, SymType, bool);
    void loadFromFileDense(std::istream&, std::vector<complex_t>&, SymType, bool);
    void loadFromFileCoo(std::istream&, std::vector<real_t>&, SymType, bool);
    void loadFromFileCoo(std::istream&, std::vector<complex_t>&, SymType, bool);

    // product of any part of dense matrix and vector common to different dense storages (child classes)
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void diagonalMatrixVector(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const; //!< diag matrix * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void diagonalVectorMatrix(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const; //!< vector * diag matrix
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void lowerMatrixVector(MatIterator&, VecIterator&, VecIterator&,
                           ResIterator&, ResIterator&, SymType = _noSymmetry) const;         //!< Lower part * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void lowerVectorMatrix(MatIterator&, VecIterator&, VecIterator&,
                           ResIterator&, ResIterator&, SymType = _noSymmetry) const;         //!< vector * lower part
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void upperMatrixVector(MatIterator&, VecIterator&, VecIterator&,
                           ResIterator&, ResIterator&, SymType = _noSymmetry) const;         //!< upper part * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void upperVectorMatrix(MatIterator&, VecIterator&, VecIterator&,
                           ResIterator&, ResIterator&, SymType = _noSymmetry) const;         //!< vector * upper part
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void rowMatrixVector(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const;   //!< row matrix * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void rowVectorMatrix(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const;   //!< vector * row matrix
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void columnMatrixVector(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const; //!< column matrix * vector
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void columnVectorMatrix(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const; //!< vector * column matrix

    //omp versions
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void parallelRowMatrixVector(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const;  //!< templated parallel row Matrix x Vector,
    template<typename MatIterator, typename V, typename R>
    void parallelRowMatrixVector(MatIterator&, const std::vector<V>&, std::vector<R>&)  const;     //!< templated parallel row Vector x Matrix, assuming square matrix
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void parallelColumnVectorMatrix(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const; //!< templated parallel column Vector x Matrix
    template<typename MatIterator, typename V, typename R>
    void parallelLowerMatrixVector(MatrixPart, MatIterator&, const std::vector<V>&, std::vector<R>&,
                                   SymType = _noSymmetry) const;                                   //!< templated parallel lower  Matrix x Vector,
    template<typename MatIterator, typename V, typename R>
    void parallelLowerVectorMatrix(MatrixPart, MatIterator&, const std::vector<V>&, std::vector<R>&,
                                   SymType = _noSymmetry) const;                                   //!< templated parallel lower Vector x Matrix,
    template<typename MatIterator, typename V, typename R>
    void parallelUpperMatrixVector(MatrixPart, MatIterator&, const std::vector<V>&, std::vector<R>&,
                                   SymType =_noSymmetry) const;                                    //!< templated parallel upper Matrix x Vector
    template<typename MatIterator, typename V, typename R>
    void parallelUpperVectorMatrix(MatrixPart, MatIterator&, const std::vector<V>&, std::vector<R>&,
                                   SymType =_noSymmetry) const;                                    //!< templated parallel upper Vector x Matrix

    void extractThreadIndex(MatrixPart, number_t&, std::vector<number_t>&) const;  //!< tool to split matrix rows/cols in balanced parts (multi-thread computation)

    //! sum of any part of dense matrix and dense matrix
    template<typename Mat1Iterator, typename Mat2Iterator, typename ResIterator>
    void sumMatrixMatrix(Mat1Iterator& itm1, Mat2Iterator& itm2, ResIterator& itrb, ResIterator& itre) const;

    //SOR and SSOR stuff
    template<typename MatIterator, typename VecIterator, typename ResIterator>
    void bzSorDiagonalMatrixVector(MatIterator&, VecIterator&, ResIterator&,
                                   ResIterator&, const real_t) const;         //!< SOR diagonal product w * D * v
    template<typename MatIterator, typename VecIterator, typename XIterator>
    void bzSorDiagonalSolver(MatIterator&, VecIterator&, XIterator&,
                             XIterator&, const real_t) const;                 //!< SOR diagonal solver D/w x = v
    template<typename MatIterator, typename VecIterator, typename XIterator>
    void bzSorLowerSolver(const MatIterator&, const MatIterator&, VecIterator&, XIterator&,
                          XIterator&, const real_t) const;                   //!< SOR lower solver (D/w+L) x = b
    template<typename MatRevIterator, typename VecRevIterator, typename XRevIterator>
    void bzSorUpperSolver(const MatRevIterator&, const MatRevIterator&, VecRevIterator&,
                          XRevIterator&, XRevIterator&, const real_t) const;  //!< SOR upper solver (D/w+U) x = b

  public:
    //Generic Gauss solver and LU factorization (see childs for faster implementation)
    template<typename T>
    void gaussSolverG(std::vector<T>&, std::vector<std::vector<T> >&) const; //!< generic Gauss solver with partial pivoting strategy  A*x1=b1, A*x2=b2, ...
    template<typename T>
    void gaussSolverG(std::vector<T>&, std::vector<T>&) const;               //!< generic Gauss solver with partial pivoting strategy  A*x=b
    template<typename T>
    void luG(std::vector<T>& A, std::vector<T>& LU, std::vector<number_t>& p) const;//!< LU factorization with row pivoting (omp)
    template<typename T>
    void luG(std::vector<T>& A, std::vector<T>& LU) const;                          //!< LU factorization with no pivoting (omp)
    //@{
    //! specializations of gaussSolver and LU factorization thats call generic functions
    virtual void gaussSolver(std::vector<real_t>& A, std::vector<std::vector<real_t> >& bs) const;
    virtual void gaussSolver(std::vector<complex_t>& A, std::vector<std::vector<complex_t> >& bs) const;
    virtual void gaussSolver(std::vector<real_t>& A, std::vector<real_t>& b) const;
    virtual void gaussSolver(std::vector<complex_t>& A, std::vector<complex_t>& b) const;
    void lu(std::vector<real_t>& A, std::vector<real_t>& LU, std::vector<number_t>& p, const SymType sym = _noSymmetry) const;
    void lu(std::vector<complex_t>& A, std::vector<complex_t>& LU, std::vector<number_t>& p, const SymType sym = _noSymmetry) const;
    void lu(std::vector<real_t>& A, std::vector<real_t>& LU, const SymType sym = _noSymmetry) const;
    void lu(std::vector<complex_t>& A, std::vector<complex_t>& LU, const SymType sym = _noSymmetry) const;
    //@}

    // Template diagonal and triangular part solvers (after factorization), call generic triangular part solvers
    //! Generic lower triangular with unit diagonal linear system solver : (I+L) x = b
    template<typename M, typename V, typename X>
    void lowerD1SolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) const;
    //@{
    //! Specializations of lower triangular part with unit diagonal linear solvers (I + L) x = v */
    void lowerD1Solver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x) const
    { lowerD1SolverG<>(m, v, x); }
    void lowerD1Solver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { lowerD1SolverG<>(m, v, x); }
    void lowerD1Solver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x) const
    { lowerD1SolverG<>(m, v, x); }
    void lowerD1Solver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { lowerD1SolverG<>(m, v, x); }
    //@}

    //! Generic lower triangular with unit diagonal linear system left solver : x (I+L) = b
    template<typename M, typename V, typename X>
    void lowerD1LeftSolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) const;
    //@{
    //! Specializations of lower triangular part with unit diagonal linear solvers (I + L) x = v */
    void lowerD1LeftSolver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x) const
    { lowerD1LeftSolverG<>(m, v, x); }
    void lowerD1LeftSolver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { lowerD1LeftSolverG<>(m, v, x); }
    void lowerD1LeftSolver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x) const
    { lowerD1LeftSolverG<>(m, v, x); }
    void lowerD1LeftSolver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { lowerD1LeftSolverG<>(m, v, x); }
    //@}

    //! Generic diagonal linear system solver : D x = b
    template<typename M, typename V, typename X>
    void diagonalSolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) const;
    //@{
    //! Specializations of diagonal linear solvers D x = v
    void diagonalSolver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x) const
    { diagonalSolverG<>(m, v, x); }
    void diagonalSolver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { diagonalSolverG<>(m, v, x); }
    void diagonalSolver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x) const
    { diagonalSolverG<>(m, v, x); }
    void diagonalSolver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { diagonalSolverG<>(m, v, x); }
    //@}

    //! Generic upper triangular with unit diagonal linear system solver : (I+U) x = b
    template<typename M, typename V, typename X>
    void upperD1SolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym = _noSymmetry) const;
    //@{
    //! Specializations of upper triangular part with unit diagonal linear solvers (I + U) x = v
    void upperD1Solver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x, const SymType sym = _noSymmetry) const
    { upperD1SolverG<>(m, v, x, sym); }
    void upperD1Solver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperD1SolverG<>(m, v, x, sym); }
    void upperD1Solver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperD1SolverG<>(m, v, x, sym); }
    void upperD1Solver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperD1SolverG<>(m, v, x, sym); }
    //@}

    //! Generic upper triangular linear system solver : (D+U) x = b
    template<typename M, typename V, typename X>
    void upperSolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym = _noSymmetry) const;
    //@{
    //! specializations of upper triangular part linear solvers (D + U) x = v
    void upperSolver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x, const SymType sym = _noSymmetry) const
    { upperSolverG<>(m, v, x, sym); }
    void upperSolver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperSolverG<>(m, v, x, sym); }
    void upperSolver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperSolverG<>(m, v, x, sym); }
    void upperSolver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperSolverG<>(m, v, x, sym); }
    //@}

    //! Generic upper triangular linear system left solver : x U = b
    template<typename M, typename V, typename X>
    void upperLeftSolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym = _noSymmetry) const;
    //@{
    //! Specializations of upper triangular linear system left solver : x U = b
    void upperLeftSolver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x, const SymType sym = _noSymmetry) const
    { upperLeftSolverG<>(m, v, x, sym); }
    void upperLeftSolver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperLeftSolverG<>(m, v, x, sym); }
    void upperLeftSolver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperLeftSolverG<>(m, v, x, sym); }
    void upperLeftSolver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const
    { upperLeftSolverG<>(m, v, x, sym); }
    //@}
};
/* --------------------------------------------------------------------------------
   Template functions
  --------------------------------------------------------------------------------*/
// print matrix of scalars
template<typename Iterator>
void DenseStorage::printScalarEntries(Iterator& it,
                                      number_t nrows, number_t ncols, number_t perRow, number_t width, number_t prec,
                                      const string_t& rowOrcol, number_t vb, std::ostream& os) const
{
  os.setf(std::ios::scientific);
  os <<eol<< "(" << words("firste") << "s " << std::min(vb, nrows) << " " << words(rowOrcol) << "s.)";
  for(number_t r = 0; r < std::min(vb, nrows) ; r++)
    {
      os << eol << "   " << rowOrcol << "   " << r + 1;
      printRowWise(os, "   ", perRow, width, prec, it, it + ncols);
      it += ncols;
    }
  os.unsetf(std::ios::scientific);
  os << std::endl;
}

template<typename Iterator>
void DenseStorage::printScalarEntriesTriangularPart(Iterator& itd, Iterator& it,
    number_t nrows, number_t ncols, number_t perRow, number_t width, number_t prec,
    const string_t& rowOrcol, number_t vb, std::ostream& os) const
{
  // itd is an iterator on matrix diagonal entries
  // it is an iterator on matrix lower or upper triangular part entries
  os.setf(std::ios::scientific);
  os << " " << words("firste") << "s " << std::min(vb, nrows) << " " << words(rowOrcol) << "s.)";
  os << eol << "   " << rowOrcol << "   " << 1 ;
  os << "   " << std::setw(width) << std::setprecision(prec) << *itd++;
  for(number_t r = 1; r < std::min(vb, nrows) ; r++)
    {
      os << eol << "   " << rowOrcol << "   " << r + 1 ;
      if(r < ncols)
        {
          printRowWise(os, "   ", perRow, width, prec, it, it + r);
          if(r % perRow == 0) { os << " ..." << eol << "..."; }
          os << std::setw(width) << std::setprecision(prec) << *itd++;
          it += r;
        }
      else
        {
          printRowWise(os, "   ", perRow, width, prec, it, it + ncols);
          it += ncols;
        }
    }
  os.unsetf(std::ios::scientific);
  os << std::endl;
}

// print matrix of matrices
template<typename Iterator>
void DenseStorage::printMatrixEntries(Iterator& it,
                                      number_t nrows, number_t ncols,
                                      const string_t& rowOrcol, number_t vb, std::ostream& os) const
{
  os.setf(std::ios::scientific);
  os << eol<<"(" << words("firste") << "s " << std::min(vb, nrows) << " " << words(rowOrcol) << "s.)";
  for(number_t r = 0; r < std::min(vb, nrows) ; r++)
    {
      os << eol << "   " << rowOrcol << "   " << r + 1;
      for(Iterator itj = it; itj != it + ncols; itj++) { os << *itj; }
      it += ncols;
    }
  os.unsetf(std::ios::scientific);
  os << std::endl;
}

template<typename Iterator>
void DenseStorage::printMatrixEntriesTriangularPart(Iterator& itd, Iterator& it,
    number_t nrows, number_t ncols,
    const string_t& rowOrcol, number_t vb, std::ostream& os) const
{
  // itd is an iterator on matrix diagonal entries
  // it is an iterator on matrix lower or upper triangular part entries
  os.setf(std::ios::scientific);
  os << " " << words("firste") << "s " << std::min(vb, nrows) << " " << words(rowOrcol) << "s.)";
  os << eol << "   " << rowOrcol << "   " << 1 ;
  os << *itd++;
  for(number_t r = 1; r < std::min(vb, nrows) ; r++)
    {
      os << eol << "   " << rowOrcol << "   " << r + 1 ;
      if(r < ncols)
        {
          for(Iterator itj = it; itj != it + r; itj++) { os << *itj; }
          os << *itd++;
          it += r;
        }
      else
        {
          for(Iterator itj = it; itj != it + ncols; itj++) { os << *itj; }
          it += ncols;
        }
    }
  os.unsetf(std::ios::scientific);
  os << std::endl;
}

/*--------------------------------------------------------------------------------
   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)
   Result vector res *** NEED NOT BE INITIALIZED HERE ***
--------------------------------------------------------------------------------*/
// partial multiplication by diagonal M*v, matrix need not be square
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::diagonalMatrixVector(MatIterator& itd, VecIterator& itvb, ResIterator& itrb, ResIterator& itre) const
{
  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 = ?!
}

// partial multiplication by diagonal v*M, matrix need not be square
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::diagonalVectorMatrix(MatIterator& itd, VecIterator& itvb, ResIterator& itrb, ResIterator& itre) const
{
  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 = ?!
}

// partial multiplication by lower triangular part
// also used in multiplication of upper part of tranposed matrix by vector
// should be called after Diagonal multiplication
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::lowerMatrixVector(MatIterator& itm, VecIterator& itvb, VecIterator& itve,
                                     ResIterator& itrb, ResIterator& itre, SymType sym) const
{
  number_t row = 1;
  number_t rox = itve - itvb;
  switch(sym)
    {
      case _skewSymmetric:
        for(ResIterator itr = itrb + 1; itr != itre; ++itr, ++row)
          for(VecIterator itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itm) { *itr -= *itm** itv; }
        break;
      case _selfAdjoint:
        for(ResIterator itr = itrb + 1; itr != itre; ++itr, ++row)
          for(VecIterator itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itm) { *itr += conj(*itm)** itv; }
        break;
      case _skewAdjoint:
        for(ResIterator itr = itrb + 1; itr != itre; ++itr, ++row)
          for(VecIterator itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itm) { *itr -= conj(*itm)** itv; }
        break;
      default:
        for(ResIterator itr = itrb + 1; itr != itre; ++itr, ++row)
          for(VecIterator itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itm) { *itr += *itm** itv; }
    }
}

// templated parallel lower Vector x Matrix, assuming square matrix
//  itm : iterator on the first element of the lower triangular part of the matrix, row access
template<typename MatIterator, typename V, typename R>
void DenseStorage::parallelLowerMatrixVector(MatrixPart lowup, MatIterator& itm, const std::vector<V>& v, std::vector<R>& r, SymType sym) const
{
#ifndef XLIFEPP_WITH_OMP
  where("DenseStorage::parallelLowerMatrixVector(...)");
  error("xlifepp_without_omp");
#else
  number_t nt = 1;
  #pragma omp parallel for lastprivate(nt)
  for(number_t i = 0; i < 1; i++) {nt = omp_get_num_threads();}

  if(nt==1) //do not use parallel computation
    {
      typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end();
      typename std::vector<R>::iterator it_rb = r.begin(), it_re = r.end();
      lowerMatrixVector(itm, itvb, itve, it_rb, it_re, sym);
      return;
    }

  number_t rox = v.size();
  typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end(), itv;
  typename std::vector<R>::iterator itrb = r.begin()+1;
  R rz0= 0.* *r.begin(), rz;
  MatIterator itmr;
  switch(sym)
    {
      case _skewSymmetric:
        #pragma omp parallel for firstprivate(itmr, itv, rz) schedule(static)
        for(number_t row=1; row<r.size(); ++row)
          {
            itmr= itm + (row*(row-1))/2;
            rz= rz0;
            for(itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itmr) { rz -= *itmr** itv;}
            *(r.begin()+ row)+=rz;
          }
        break;
      case _selfAdjoint:
        #pragma omp parallel for firstprivate(itmr, itv, rz) schedule(static)
        for(number_t row=1; row<r.size(); ++row)
          {
            itmr= itm + (row*(row-1))/2;
            rz= rz0;
            for(itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itmr) { rz += conj(*itmr)** itv; }
            *(r.begin()+ row)+=rz;
          }
        break;
      case _skewAdjoint:
        #pragma omp parallel for firstprivate(itmr, itv, rz) schedule(static)
        for(number_t row=1; row<r.size(); ++row)
          {
            itmr= itm + (row*(row-1))/2;
            rz= rz0;
            for(itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itmr) { rz -= conj(*itmr)** itv; }
            *(r.begin()+ row)+=rz;
          }
        break;
      default:
        #pragma omp parallel for firstprivate(itmr, itv, rz) schedule(dynamic)
        for(number_t row=1; row<r.size(); ++row)
          {
            itmr= itm + (row*(row-1))/2;
            rz= rz0;
            for(itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itmr) { rz += *itmr** itv; }
            *(r.begin()+ row)+=rz;
          }
    }
#endif
}

// not included parallel case, see parallelLowerVectorMatrix for its parallelized version
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::lowerVectorMatrix(MatIterator& itm, VecIterator& itvb, VecIterator& itve,
                                     ResIterator& itrb, ResIterator& itre, SymType sym) const
{
  number_t col = 1;
  number_t cox = itre - itrb;
  switch(sym)
    {
      case _skewSymmetric:
        for(VecIterator itv = itvb + 1; itv != itve; ++itv, ++col)
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) { *itr -= *itv** itm ; }
        break;
      case _selfAdjoint:
        for(VecIterator itv = itvb + 1; itv != itve; ++itv, ++col)
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) { *itr += *itv * conj(*itm); }
        break;
      case _skewAdjoint:
        for(VecIterator itv = itvb + 1; itv != itve; ++itv, ++col)
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) { *itr -= *itv * conj(*itm); }
        break;
      default:
        for(VecIterator itv = itvb + 1; itv != itve ; itv++, col++)
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) *itr += *itv** itm;
    }
}

// templated parallel lower Vector x Matrix, assuming square matrix
//  itm : iterator on the first element of the lower triangular part of the matrix, row access
template<typename MatIterator, typename V, typename R>
void DenseStorage::parallelLowerVectorMatrix(MatrixPart lowup, MatIterator& itm, const std::vector<V>& v, std::vector<R>& r, SymType sym) const
{
#ifndef XLIFEPP_WITH_OMP
  where("DenseStorage::parallelLowerVectorMatrix(...)");
  error("xlifepp_without_omp");
#else
  number_t nt = 1;
  std::vector<number_t> threadIndex;
  extractThreadIndex(lowup,nt,threadIndex);
  if(nt==1) //do not use parallel computation
    {
      typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end();
      typename std::vector<R>::iterator it_rb = r.begin(), it_re = r.end();
      lowerVectorMatrix(itm, itvb, itve, it_rb, it_re, sym);
      return;
    }
  //compute product by distributing matrix parts given by threadIndex vector on different thread
  std::vector<std::vector<R> > res(nt,std::vector<R>(r.size(),0.*r[0])); //temporary vector to avoid data races
  number_t cox=r.size();
  #pragma omp parallel for
  for(number_t t=0; t<nt; ++t)  //partial product in parallel
    {
      number_t r1 = threadIndex[t], r2 = threadIndex[t+1];
      MatIterator itmt = itm + (r1*(r1+1)/2);
      typename std::vector<V>::const_iterator itvb = v.begin()+ (r1+1), itve = v.begin() + (r2+1), itv;
      if(t==nt-1) itve=v.end();
      typename std::vector<R>::iterator itrb = res[t].begin(), itre = res[t].end(), itr;
      number_t col = r1+1;
      switch(sym)
        {
          case _skewSymmetric:
            for(itv = itvb; itv != itve; ++itv, ++col)
              for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr -= *itv** itmt ; }
            break;
          case _selfAdjoint:
            for(itv = itvb; itv != itve; ++itv, ++col)
              for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr += *itv * conj(*itmt); }
            break;
          case _skewAdjoint:
            for(itv = itvb; itv != itve; ++itv, ++col)
              for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr -= *itv * conj(*itmt); }
            break;
          default:
            for(itv = itvb; itv != itve ; itv++, col++)
              for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr += *itv** itmt; }
        }
    }
  //accumulate temporary result in a serial way
  typename std::vector<R>::iterator itrb = r.begin(), itre = r.end(), itr;
  typename std::vector<R>::iterator itrp;
  for(number_t t=0; t<nt; ++t)
    {
      itrp = res[t].begin();
      for(itr=itrb; itr!=itre; ++itr, ++itrp) *itr+=*itrp;
    }
#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
// see parallelUpperMatrixVector for its parallelized version
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::upperMatrixVector(MatIterator& itm, VecIterator& itvb, VecIterator& itve,
                                     ResIterator& itrb, ResIterator& itre, SymType sym) const
{
  number_t col = 1;
  number_t cox = itre - itrb;
  switch(sym)
    {
      case _skewSymmetric:
        for(VecIterator itv = itvb + 1; itv != itve; ++itv, ++col)
          //for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); itr++, itm++) { *itr -= *itm * *itv; }
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) { *itr -= *itv** itm ; }
        break;
      case _selfAdjoint:
        for(VecIterator itv = itvb + 1; itv != itve; itv++, col++)
          //for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); itr++, itm++) { *itr += conj(*itm) * *itv; }
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) { *itr += *itv * conj(*itm); }
        break;
      case _skewAdjoint:
        for(VecIterator itv = itvb + 1; itv != itve; itv++, ++col)
          //for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); itr++, itm++) { *itr -= conj(*itm) * *itv; }
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) { *itr -= *itv * conj(*itm); }
        break;
      case _symmetric:
        for(VecIterator itv = itvb + 1; itv != itve; itv++, ++col)
          //for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); itr++, itm++) { *itr += *itm * *itv; }
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) { *itr += *itv** itm; }
        break;
      default:
        for(VecIterator itv = itvb + 1; itv != itve; itv++, col++)
          for(ResIterator itr = itrb;  itr != itrb + std::min(col, cox); ++itr, ++itm) { *itr += *itm** itv; }
    }
}

//templated parallel upper Matrix x Vector, assuming square matrix
// itm : iterator on the first element of the upper triangular part of the matrix, column access
template<typename MatIterator, typename V, typename R>
void DenseStorage::parallelUpperMatrixVector(MatrixPart lowup, MatIterator& itm, const std::vector<V>& v, std::vector<R>& r, SymType sym) const
{
#ifdef XLIFEPP_WITH_OMP
  number_t nt = 1;
  std::vector<number_t> threadIndex;
  extractThreadIndex(lowup,nt,threadIndex);
  if(nt==1) //do not use parallel computation
    {
      typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end();
      typename std::vector<R>::iterator it_rb = r.begin(), it_re = r.end();
      upperMatrixVector(itm, itvb, itve, it_rb, it_re, sym);
      return;
    }
  //compute product by distributing matrix parts given by threadIndex vector on different thread
  std::vector<std::vector<R> > res(nt,std::vector<R>(r.size(),0.*r[0])); //temporary vector to avoid data races
  number_t cox=r.size();
  #pragma omp parallel for
  for(number_t t=0; t<nt; ++t)  //partial product in parallel
    {
      number_t r1 = threadIndex[t], r2 = threadIndex[t+1];
      MatIterator itmt = itm + (r1*(r1+1)/2);
      typename std::vector<V>::const_iterator itvb = v.begin()+ (r1+1), itve = v.begin() + (r2+1), itv;
      if(t==nt-1) itve=v.end();
      typename std::vector<R>::iterator itrb = res[t].begin(), itre = res[t].end(), itr;
      number_t col = r1+1;
      switch(sym)
        {
          case _skewSymmetric:
            for(itv = itvb; itv != itve; ++itv, ++col)
              //for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr -= *itv** itmt ; }
              for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr -= *itmt * *itv; }
            break;
          case _selfAdjoint:
            for(itv = itvb; itv != itve; ++itv, ++col)
              //for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr += *itv * conj(*itmt); }
              for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr += conj(*itmt) * *itv; }
            break;
          case _skewAdjoint:
            for(itv = itvb; itv != itve; ++itv, ++col)
              //for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr -= *itv * conj(*itmt); }
              for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr -=  conj(*itmt) ** itv; }
            break;
          default:
            for(itv = itvb; itv != itve ; itv++, col++)
              //for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr += *itv** itmt; }
              for(itr = itrb;  itr != itrb + std::min(col,cox); ++itr, ++itmt) { *itr += * itmt * *itv; }
        }
    }
  //accumulate temporary result in a serial way
  typename std::vector<R>::iterator itrb = r.begin(), itre = r.end(), itr;
  typename std::vector<R>::iterator itrp;
  for(number_t t=0; t<nt; ++t)
    {
      itrp = res[t].begin();
      for(itr=itrb; itr!=itre; ++itr, ++itrp) *itr+=*itrp;
    }
#endif
}

template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::upperVectorMatrix(MatIterator& itm, VecIterator& itvb, VecIterator& itve,
                                     ResIterator& itrb, ResIterator& itre, SymType sym) const
{
  number_t col= 1;
  number_t cox = itve - itvb;
  switch(sym)
    {
      case _skewSymmetric:
        for(ResIterator itr = itrb + 1; itr != itre; ++itr, ++col)
          for(VecIterator itv = itvb; itv != itvb + std::min(col, cox); ++itv, ++itm) { *itr -= *itm** itv; }
        //for(VecIterator itv = itvb; itv != itvb + std::min(col, cox); ++itv, ++itm) { *itr -= *itv **itm;} //vector of vector
        break;
      case _selfAdjoint:
        for(ResIterator itr = itrb + 1; itr != itre; ++itr, ++col)
          for(VecIterator itv = itvb; itv != itvb + std::min(col, cox); ++itv, ++itm) { *itr += conj(*itm)** itv; }
        //for(VecIterator itv = itvb; itv != itvb + std::min(col, cox); ++itv, ++itm) { *itr += *itv * conj(*itm);} //vector of vector
        break;
      case _skewAdjoint:
        for(ResIterator itr = itrb + 1; itr != itre; ++itr, ++col)
          for(VecIterator itv = itvb; itv != itvb + std::min(col, cox); ++itv, ++itm) { *itr -= conj(*itm)** itv; }
        //for(VecIterator itv = itvb; itv != itvb + std::min(col, cox); ++itv, ++itm) { *itr -= *itv * conj(*itm);} //vector of vector
        break;
      default:
        for(ResIterator itr = itrb + 1; itr != itre; ++itr, ++col)
          for(VecIterator itv = itvb; itv != itvb + std::min(col, cox); ++itv, ++itm) { *itr += *itm** itv; }
        //for(VecIterator itv = itvb; itv != itvb + std::min(col, cox); ++itv, ++itm) { *itr += *itv **itm;} //vector of vector
    }
}

// templated parallel lower Vector x Matrix, assuming square matrix
//  itm : iterator on the first element of the lower triangular part of the matrix, row access
template<typename MatIterator, typename V, typename R>
void DenseStorage::parallelUpperVectorMatrix(MatrixPart lowup, MatIterator& itm, const std::vector<V>& v, std::vector<R>& r, SymType sym) const
{
#ifndef XLIFEPP_WITH_OMP
  where("DenseStorage::parallelLowerMatrixVector(...)");
  error("xlifepp_without_omp");
#else
  number_t nt = 1;
  std::vector<number_t> threadIndex;
  extractThreadIndex(lowup,nt,threadIndex);
  if(nt==1) //do not use parallel computation
    {
      typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end();
      typename std::vector<R>::iterator it_rb = r.begin(), it_re = r.end();
      lowerMatrixVector(itm, itvb, itve, it_rb, it_re, sym);
      return;
    }

  //multi-threads
  number_t rox = v.size();
  typename std::vector<V>::const_iterator itvb = v.begin(), itve = v.end(), itv;
  typename std::vector<R>::iterator itrb = r.begin()+1;
  R rz0 = 0.* *r.begin();
  switch(sym)
    {
      case _skewSymmetric:
        #pragma omp parallel for schedule(static)
        for(number_t row=1; row<r.size(); ++row)
          {
            typename std::vector<V>::const_iterator itv;
            MatIterator itmr= itm + (row*(row-1)/2);
            R rz= rz0;
            for(itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itmr) { rz -= * itv * *itmr; }
            typename std::vector<R>::iterator itr = r.begin()+ row;
            *itr+=rz;
          }
        break;
      case _selfAdjoint:
        #pragma omp parallel for schedule(static)
        for(number_t row=1; row<r.size(); ++row)
          {
            typename std::vector<V>::const_iterator itv;
            MatIterator itmr= itm + (row*(row-1)/2);
            R rz= rz0;
            for(itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itmr) { rz += * itv * conj(*itmr); }
            typename std::vector<R>::iterator itr = r.begin()+ row;
            *itr+=rz;
          }
        break;
      case _skewAdjoint:
        #pragma omp parallel for schedule(static)
        for(number_t row=1; row<r.size(); ++row)
          {
            typename std::vector<V>::const_iterator itv;
            MatIterator itmr= itm + (row*(row-1)/2);
            R rz= rz0;
            for(itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itmr) { rz -= * itv * conj(*itmr); }
            typename std::vector<R>::iterator itr = r.begin()+ row;
            *itr+=rz;
          }
        break;
      default:
        #pragma omp parallel for schedule(static)
        for(number_t row=1; row<r.size(); ++row)
          {
            typename std::vector<V>::const_iterator itv;
            MatIterator itmr= itm + (row*(row-1)/2);
            R rz= rz0;
            for(itv = itvb; itv != itvb + std::min(row, rox); ++itv, ++itmr) { rz += * itv * *itmr; }
            typename std::vector<R>::iterator itr = r.begin()+ row;
            *itr+=rz;
          }
    }
#endif
}

// multiplication M x V of dense 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 DenseStorage::rowMatrixVector(MatIterator& itm, VecIterator& itvb, VecIterator& itve, ResIterator& itrb, ResIterator& itre)  const
{
  for(ResIterator itr = itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  for(ResIterator itr = itrb; itr != itre; ++itr)
    for(VecIterator itv = itvb; itv != itve; ++itv, ++itm) { *itr += *itm** itv; }
}

// also used in multiplication of column-major transposed matrix by vector
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::parallelRowMatrixVector(MatIterator& itm, VecIterator& itvb, VecIterator& itve, ResIterator& itrb, ResIterator& itre)  const
{
  for(ResIterator itr = itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
#ifndef XLIFEPP_WITH_OMP
  for(ResIterator itr = itrb; itr != itre; ++itr)
    for(VecIterator itv = itvb; itv != itve; ++itv, ++itm) { *itr += *itm** itv; }
#else
  number_t nr =nbOfRows(), nc=nbOfColumns();
  #pragma omp parallel for
  for(number_t r=0; r<nr; ++r)
    {
      ResIterator itr = itrb + r;
      MatIterator itmr =  itm + r*nc;
      for(VecIterator itv = itvb; itv != itve; ++itv, ++itmr) { *itr += *itmr** itv; }
    }
#endif // XLIFEPP_WITH_OMP
}

// also used in multiplication of column-major transposed matrix by vector
template<typename MatIterator, typename V, typename R>
void DenseStorage::parallelRowMatrixVector(MatIterator& itm, const std::vector<V>& v, std::vector<R>& r)  const
{
#ifndef XLIFEPP_WITH_OMP
  where("DenseStorage::parallelRowMatrixVector(...)");
  error("xlifepp_without_omp");
#else
  typename std::vector<R>::iterator itrb = r.begin(), itre=r.end(), itr;
  for(itr=itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  number_t nr =nbOfRows(), nc=nbOfColumns();
  typename std::vector<V>::const_iterator itvb = v.begin(), itve=v.end(), itv;
  MatIterator itmr ;
  R t;
  #pragma omp parallel for firstprivate(itr, itv, itmr, t) schedule(static)
  for(number_t r=0; r<nr; ++r)
    {
      itr = itrb + r;
      itmr =  itm + r*nc;
      t = *itr;
      for(itv = itvb; itv != itve; ++itv, ++itmr) { t += *itmr** itv; }
      *itr = t;
    }
#endif // XLIFEPP_WITH_OMP
}

//template<typename MatIterator, typename V, typename R>
//void DenseStorage::parallelRowMatrixVector(MatIterator& itm, const std::vector<V>& v, std::vector<R>& r)  const
//{
//#ifndef XLIFEPP_WITH_OMP
//  where("DenseStorage::parallelRowMatrixVector(...)");
//  error("xlifepp_without_omp");
//#else
//  typename std::vector<R>::iterator itrb = r.begin(), itre=r.end(), itr;
//  for(itr=itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
//  number_t nr =nbOfRows(), nc=nbOfColumns();
//  typename std::vector<V>::const_iterator itvb = v.begin(), itve=v.end(), itv;
//  MatIterator itmr ;
//  R t;
//  number_t br=0, rm, sr=100;  //size of block
//  while(br<nr)
//  {
//   rm=std::min(br+sr,nr);
//   #pragma omp parallel for firstprivate(itr, itv, itmr, t) schedule(static)
//   for(number_t r=br; r < rm; ++r)
//    {
//      itr = itrb + r;
//      itmr =  itm + r*nc;
//      t = *itr;
//      for(itv = itvb; itv != itve; ++itv, ++itmr) { t += *itmr** itv; }
//      *itr = t;
//    }
//    br+=sr;
//  }
// #endif // XLIFEPP_WITH_OMP
//}

//template<typename MatIterator, typename V, typename R>
//void DenseStorage::parallelRowMatrixVector(MatIterator& itm, const std::vector<V>& v, std::vector<R>& r)  const
//{
//#ifndef XLIFEPP_WITH_OMP
//  where("DenseStorage::parallelRowMatrixVector(...)");
//  error("xlifepp_without_omp");
//#else
//  typename std::vector<R>::iterator itrb = r.begin(), itre=r.end(), itr;
//  for(itr=itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
//  number_t nr =nbOfRows(), nc=nbOfColumns();
//  typename std::vector<V>::const_iterator itvb = v.begin(), itve=v.end(), itv;
//  MatIterator itmr ;
//  R t;
//  number_t sr=100;  //size of block
//  number_t nbr=nr/sr, rm;
//  #pragma omp parallel for firstprivate(itr, itv, itmr, t) schedule(static)
//  for(number_t br=0; br<=nbr; ++br)
//  {
//   number_t r0=br*sr, rm=std::min(r0+sr,nr);
//   for(number_t r=r0; r < rm; ++r)
//    {
//      itr = itrb + r;
//      itmr =  itm + r*nc;
//      t = *itr;
//      for(itv = itvb; itv != itve; ++itv, ++itmr) { t += *itmr** itv; }
//      *itr = t;
//    }
//  }
// #endif // XLIFEPP_WITH_OMP
//}

// multiplication V x M of vector V and dense row-major access matrix M
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::rowVectorMatrix(MatIterator& itm, VecIterator& itvb, VecIterator& itve, ResIterator& itrb, ResIterator& itre)  const
{
  for(ResIterator itr = itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  for(VecIterator itv = itvb; itv != itve; ++itv)
    for(ResIterator itr = itrb; itr != itre; ++itr, ++itm) { *itr += *itv** itm; }
}

// multiplication M x V of dense 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 DenseStorage::columnMatrixVector(MatIterator& itm, VecIterator& itvb, VecIterator& itve, ResIterator& itrb, ResIterator& itre) const
{
  for(ResIterator itr = itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  for(VecIterator itv = itvb; itv != itve; ++itv)
    for(ResIterator itr = itrb; itr != itre; ++itr, ++itm) { *itr += *itm** itv; }
}

// multiplication V x M of vector V and dense column-major access matrix M
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::columnVectorMatrix(MatIterator& itm, VecIterator& itvb, VecIterator& itve, ResIterator& itrb, ResIterator& itre) const
{
  for(ResIterator itr = itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
  for(ResIterator itr = itrb; itr != itre; ++itr)
    for(VecIterator itv = itvb; itv != itve; ++itv, ++itm) *itr += *itv** itm;
  // Transpose (*itm) HERE ???????????????????????
}

// multiplication V x M of vector V and dense column-major access matrix M
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::parallelColumnVectorMatrix(MatIterator& itm, VecIterator& itvb, VecIterator& itve, ResIterator& itrb, ResIterator& itre) const
{
  for(ResIterator itr = itrb; itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
#ifndef XLIFEPP_WITH_OMP
  for(ResIterator itr = itrb; itr != itre; ++itr)
    for(VecIterator itv = itvb; itv != itve; ++itv, ++itm) *itr += *itv** itm;
#else
  number_t nr=nbOfRows(), nc=nbOfColumns();
  #pragma omp parallel for
  for(number_t c=0; c<nc; ++c)
    {
      ResIterator itr = itrb + c;
      MatIterator itmr =  itm + c*nr;
      for(VecIterator itv = itvb; itv != itve; ++itv, ++itmr) { *itr += *itmr** itv; }
    }
#endif
  // Transpose (*itm) HERE ???????????????????????
}


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

/* -------------------------------------------------------------------------------------------------------
                                   generic lower and upper solvers
   -------------------------------------------------------------------------------------------------------*/
// solve LD1 * x = v
template<typename M, typename V, typename X>
void DenseStorage::lowerD1SolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) const
{
    typename std::vector<V>::iterator itv=v.begin();
    typename std::vector<X>::iterator itxb=x.begin(), itx;
    number_t n=x.size();
    for(number_t k=1; k<=n; ++k, ++itv)
    {
        X t=*itv;
        itx=itxb;
        for(number_t j=1; j<k; ++j,++itx)
        {
            t -= m[pos(k,j)] * *itx;
        }
        *itx = t;
    }
}

// solve U * x = v
template<typename M, typename V, typename X>
void DenseStorage::upperSolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym) const
{
    number_t n=x.size();
    typename std::vector<V>::reverse_iterator itv=v.rbegin()++;
    typename std::vector<X>::reverse_iterator itxb=x.rbegin()++, itx;
    for(number_t k=n; k>0; k--, ++itv)
    {
        switch(sym)
        {
            case _skewSymmetric :
            {
                X t=*itv; itx=itxb;
                for(number_t j=n; j>k; j--, ++itx) t += m[pos(k, j, sym)] * *itx;
                *itx = t / m[pos(k,k)];
            }
            break;
            case _selfAdjoint :
            {
                X t=*itv; itx=itxb;
                for(number_t j=n; j>k; j--, ++itx) t -= conj(m[pos(k, j, sym)]) * *itx;
                *itx = t / m[pos(k,k)];
            }
            break;
            case _skewAdjoint :
            {
                X t=*itv; itx=itxb;
                for(number_t j=n; j>k; j--, ++itx) t += conj(m[pos(k, j, sym)]) * *itx;
                *itx = t / m[pos(k,k)];
            }
            break;
            default : //symmetric
                {
                X t=*itv; itx=itxb;
                for(number_t j=n; j>k; j--, ++itx) t -= m[pos(k, j, sym)] * *itx;
                *itx = t / m[pos(k,k)];
                }
        }
    }
}

// solve UD1 * x = v
template<typename M, typename V, typename X>
void DenseStorage::upperD1SolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym) const
{
    number_t n=x.size();
    typename std::vector<V>::reverse_iterator itv=v.rbegin()++;
    typename std::vector<X>::reverse_iterator itxb=x.rbegin()++, itx;
    for(number_t k=n; k>0; k--, ++itv)
    {
       switch(sym)
        {
            case _skewSymmetric :
            {
                X t=*itv; itx=itxb;
                for(number_t j=n; j>k; j--, ++itx) t += m[pos(k, j, sym)] * *itx;
                *itx = t;
            }
            break;
            case _selfAdjoint :
            {
                X t=*itv; itx=itxb;
                for(number_t j=n; j>k; j--, ++itx) t -= conj(m[pos(k, j, sym)]) * *itx;
                *itx = t ;
            }
            break;
            case _skewAdjoint :
            {
                X t=*itv; itx=itxb;
                for(number_t j=n; j>k; j--, ++itx) t += conj(m[pos(k, j, sym)]) * *itx;
                *itx = t;
            }
            break;
            default : //symmetric
                {
                X t=*itv; itx=itxb;
                for(number_t j=n; j>k; j--, ++itx) t -= m[pos(k, j, sym)] * *itx;
                *itx = t;
                }
        }
    }
}

// solve x * LD1 = v
template<typename M, typename V, typename X>
void DenseStorage::lowerD1LeftSolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) const
{
    number_t n=x.size();
    typename std::vector<V>::reverse_iterator itv=v.rbegin()++;
    typename std::vector<X>::reverse_iterator itxb=x.rbegin()++, itx;
    for(number_t k=n; k>0; k--, ++itv)
    {
        X t=*itv; itx=itxb;
        for(number_t j=n; j>k; j--, ++itx) t -= *itx * m[pos(j, k)] ;
        *itx = t;
    }
}

// solve x * U = v
template<typename M, typename V, typename X>
void DenseStorage::upperLeftSolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym) const
{
    typename std::vector<V>::iterator itv=v.begin();
    typename std::vector<X>::iterator itxb=x.begin(), itx;
    number_t n=x.size();
    for(number_t k=1; k<=n; ++k, ++itv)
    {
        switch(sym)
        {
            case _skewSymmetric :
            {
                X t=*itv; itx=itxb;
                for(number_t j=1; j<k; ++j,++itx) t += *itx * m[pos(j, k, sym)];
                *itx = t/m[pos(k,k)];
            }
            break;
            case _selfAdjoint :
            {
                X t=*itv; itx=itxb;
                for(number_t j=1; j<k; ++j,++itx) t -= *itx * conj(m[pos(j, k, sym)]);
                *itx = t/m[pos(k,k)];
            }
            break;
            case _skewAdjoint :
            {
                X t=*itv; itx=itxb;
                for(number_t j=1; j<k; ++j,++itx) t += *itx * conj(m[pos(j, k, sym)]);
                *itx = t/m[pos(k,k)];
            }
            break;
            default : //symmetric
                {
                X t=*itv; itx=itxb;
                for(number_t j=1; j<k; ++j,++itx) t -= *itx * m[pos(j, k, sym)];
                *itx = t/m[pos(k,k)];
                }
        }
    }
}

// solve D * x = v
template<typename M, typename V, typename X>
void DenseStorage::diagonalSolverG(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) const
{
    typename std::vector<V>::iterator itv=v.begin();
    typename std::vector<X>::iterator itx=x.begin();
    number_t n=x.size();
    for(number_t k=1; k<=n; k++, ++itv, ++itx)
    {
        *itx = *itv / m[pos(k,k)] ;
    }
}

/* -------------------------------------------------------------------------------------------------------
                                    SOR and SSOR stuff
   -------------------------------------------------------------------------------------------------------*/
/*! SOR-like diagonal product  r = w * D * v
    \param itd is a (forward) iterator on matrix diagonal entries (D)
    \param itvb is a (forward) iterator on right hand side vector v entries
    \param itrb, itre are (forward) iterators on solution vector r entries
    \param w is the SOR parameter
*/
template<typename MatIterator, typename VecIterator, typename ResIterator>
void DenseStorage::bzSorDiagonalMatrixVector(MatIterator& itd, VecIterator& itvb, ResIterator& itrb,
    ResIterator& itre, const real_t w) const
{
  VecIterator itv = itvb;
  ResIterator itr;
  for(itr = itrb; itr != itrb + diagonalSize(); ++itr, ++itd, ++itv) { *itr = w** itv** itd; }
  for(itr = itrb + diagonalSize(); itr != itre; ++itr) { *itr *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
}

/*! SOR-like diagonal D/w x = v solver
    \param itd  is a (forward) iterator on matrix diagonal entries (D)
    \param itvb is a (forward) iterator on right hand side vector v entries
    \param itxb, itxe are (forward) iterators on solution vector x entries
    \param w is the SOR parameter
*/
template<typename MatIterator, typename VecIterator, typename XIterator>
void DenseStorage::bzSorDiagonalSolver(MatIterator& itd, VecIterator& itvb, XIterator& itxb,
                                       XIterator& itxe, const real_t w) const
{
  VecIterator itv = itvb;
  XIterator itx;
  for(itx = itxb; itx != itxb + diagonalSize(); ++itx, ++itd, ++itv) { *itx = w** itv / *itd; }
  for(itx = itxb + diagonalSize(); itx != itxe; ++itx) { *itx *= 0; }     // *itr=0 : ambiguous resolution in Vector<K>::operator = ?!
}

/*! 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 w is the SOR parameter
*/
template<typename MatIterator, typename VecIterator, typename XIterator>
void DenseStorage::bzSorLowerSolver(const MatIterator& itdb, const MatIterator& itlb, VecIterator& itbb, XIterator& itxb,
                                    XIterator& itxe, const real_t w) const
{
  MatIterator itd=itdb, itl=itlb;
  XIterator itx = itxb, itxj;
  number_t row=0;
  for(itx = itxb; itx != itxe; ++itx, ++row, ++itd)
    {
      *itx = *(itbb + row);
      itxj=itxb;
      for(number_t j=0; j<row; ++j, ++itl, ++itxj)  *itx -= *itl** itxj;
      *itx *= w / *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 "++"
    \param itrdb is a (backward) iterator on matrix diagonal entries (D), starting on matrix last diagonal entry
    \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 w is the SOR parameter
*/
template<typename MatRevIterator, typename VecRevIterator, typename XRevIterator>
void DenseStorage::bzSorUpperSolver(const MatRevIterator& itrdb, const MatRevIterator& itrub, VecRevIterator& itrbb,
                                    XRevIterator& itrxb, XRevIterator& itrxe, const real_t w) const
{
  MatRevIterator itrd=itrdb, itru=itrub;
  VecRevIterator itrb = itrbb;
  XRevIterator itrx = itrxb, itrxi;
  for(itrx = itrxb; itrx != itrxe; ++itrx) { *itrx = *itrb++;}  //initialize x with b

  itrx = itrxb;
  number_t nbc=nbOfColumns();
  for(number_t c=nbc; c!=0; --c, ++itrd, ++itrx)
    {
      *itrx *= w / *itrd;
      itrxi = itrxb + (nbc -c + 1);
      for(number_t i=1; i< c; ++i, ++itrxi, ++itru) *itrxi -= *itru** itrx;
    }
}

} // end of namespace xlifepp

#endif // DENSE_STORAGE_HPP