largeMatrix/skylineStorage/DualSkylineStorage.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 DualSkylineStorage.hpp
  \authors D. Martin, E. Lunéville, NGUYEN Manh Ha, N. Kielbasiewicz
  \since 24 juin 2011
  \date 6 Juin 2014

  \brief Definition of the xlifepp::DualSkylineStorage class

  xlifepp::DualSkylineStorage class is the child class of xlifepp::SkylineStorage class dealing with skyline Row/Col Storage

  The diagonal is always first stored, then the lower triangular part in row skyline, then the upper triangular part in col skyline
  NOTE : the first entry (say (1,1)) is stored at adress 1 in values vector (not at adress 0!)

  - rowPointer_    : for lower part, vector of positions of begining of rows in previous vector
  - colPointer_    : for upper part, vector of positions of begining of columns in previous vector


  Example : The following 5 x 6 non symmetric matrix
  -------
  0  d x . x . x
  1  . d x . . x     ndiag       = 5
  2  x x d . . .     rowPointer_ = [ 0 0 0 2 3 10 ]     (last gives the number of non zeros in the strict lower part)
  3  . . x d . .     colPointer_ = [ 0 0 1 2 5 5 10 ]   (last gives the number of non zeros in the strict upper part)
  4  x . x x d .

  Note :  for an empty row, rowPointer_is equal to the previous one : length of row i = rowPointer_[i+1]-rowPointer_[i]
          for an empty col, colPointer_is equal to the previous one : length of col j = colPointer_[j+1]-colPointer_[j]
*/

#ifndef DUAL_SKYLINE_STORAGE_HPP
#define DUAL_SKYLINE_STORAGE_HPP

#include "config.h"
#include "SkylineStorage.hpp"

namespace xlifepp
{

/*!
   \class DualSkylineStorage
   child class for dual row/col compressed sparse storage
 */

class DualSkylineStorage : public SkylineStorage
{
  protected:
    std::vector<number_t> rowPointer_; //!< vector of positions of begining of rows
    std::vector<number_t> colPointer_; //!< vector of positions of begining of cols

  public:
    // constructor, destructor
    DualSkylineStorage(number_t nr = 0, number_t nc = 0, string_t id="DualSkylineStorage");  //!< default constructor
    DualSkylineStorage(const std::vector<number_t>&, const std::vector<number_t>&,
                       string_t  id="DualSkylineStorage");                               //!< explicit constructor
    ~DualSkylineStorage() {}                                                           //!< destructor
    //! constructor by a pair of global numerotation vectors
    DualSkylineStorage(number_t, number_t, const std::vector< std::vector<number_t> >&, const std::vector< std::vector<number_t> >&, string_t id="DualSkylineStorage");
    //! constructor by the list of col index by rows
    template<class L>
    DualSkylineStorage(number_t, number_t, const std::vector<L>&, string_t id="DualSkylineStorage");
    DualSkylineStorage* clone() const                                 //! create a clone (virtual copy constructor, covariant)
    {return new DualSkylineStorage(*this);}
    DualSkylineStorage* toScalar(dimen_t, dimen_t);                       //!< create a new scalar DualSkyline storage from current DualSkyline storage and submatrix sizes
    virtual void clear()                                              //! clear storage vectors
    {rowPointer_.clear();colPointer_.clear();}

    // utilities
    number_t size() const                                           //! number of non zero elements stored
    {return std::min(nbRows_, nbCols_) + rowPointer_[nbRows_] + colPointer_[nbCols_];}
    number_t lowerPartSize() const {return rowPointer_[nbRows_];}   //!< returns number of entries of lower triangular part
    number_t upperPartSize() const {return colPointer_[nbCols_];}   //!< returns number of entries of upper triangular part
    const std::vector<number_t>& rowPointer() const                 //! returns rowPointer_ (const)
    {return rowPointer_;}
    const std::vector<number_t>& colPointer() const                 //! returns colPointer_ (const)
    {return colPointer_;}
    bool sameStorage(const MatrixStorage&) const;                   //!< check if two storages have the same structures

    // access operator
    number_t pos(number_t i, number_t j, SymType s = _noSymmetry) const;     //!< returns adress of entry (i,j)
    void positions(const std::vector<number_t>&, const std::vector<number_t>&,
                   std::vector<number_t>&, bool errorOn = true, SymType = _noSymmetry) const; //!< access to submatrix positions
    std::vector<std::pair<number_t, number_t> > getCol(SymType s, number_t c, number_t r1=1, number_t r2=0) const; //!< get (row indices, adress) of col c in set [r1,r2]
    std::vector<std::pair<number_t, number_t> > getRow(SymType s, number_t r, number_t c1=1, number_t c2=0) const; //!< get (col indices, adress) of row r in set [c1,c2]
    void addSubMatrixIndices(const std::vector<number_t>&,const std::vector<number_t>&); //!< add dense submatrix indices to storage

    /* template specializations
       This rather COMPLICATED implementation of multMatrixVector is made
       to simulate "template virtual functions" in abstract base class MatrixStorage
       and thus to avoid "switch...case:" programming and "downcasting" in Matrix*Vector
       overloaded operator* of class LargeMatrix.
       Note that the call to template functions above is forced with the "<>" template
       descriptor in the following specialized virtual functions (this is not recursion) */

    /*------------------------------------------------------------------------------------------------------------------
                                        input and output stuff
      ------------------------------------------------------------------------------------------------------------------*/
    void print(std::ostream&) const;                                     //!< print storage pointers
    void printEntries(std::ostream&, const std::vector<real_t>&,
                      number_t vb, const SymType sym) const;             //!< print real scalar matrix
    void printEntries(std::ostream&, const std::vector<complex_t>&,
                      number_t vb, const SymType sym) const;             //!< print complex scalar matrix
    void printEntries(std::ostream&, const std::vector<Matrix<real_t> >&,
                      number_t vb, const SymType sym) const;             //!< print matrix of real matrices
    void printEntries(std::ostream&, const std::vector<Matrix<complex_t> >&,
                      number_t vb, const SymType sym) const;             //!< print matrix of complex matrices
    void printEntries(std::ostream&, const std::vector<Vector<real_t> >&,
                      number_t vb, const SymType sym) const;             //!< print matrix of real vectors (not available)
    void printEntries(std::ostream&, const std::vector<Vector<complex_t> >&,
                      number_t vb, const SymType sym) const;             //!< print matrix of complex vectors (not available)

    template<typename MatIterator>
    void printCooMatrix(MatIterator&, SymType, std::ostream&) const; //!< output matrix of real scalars in coordinate form
    void printCooMatrix(std::ostream&, const std::vector<real_t>&,
                        SymType s = _noSymmetry) const;       //!< output matrix of real scalars in coordinate form
    void printCooMatrix(std::ostream&, const std::vector<complex_t>&,
                        SymType s = _noSymmetry) const;       //!< output matrix of real scalars in coordinate form
    void printCooMatrix(std::ostream&, const std::vector< Matrix<real_t> >&,
                        SymType s = _noSymmetry) const;       //!< output matrix of real scalars in coordinate form
    void printCooMatrix(std::ostream&, const std::vector< Matrix<complex_t> >&,
                        SymType s = _noSymmetry) const;       //!< output matrix of complex scalars in coordinate form

    void print(PrintStream& os) const
      {print(os.currentStream());}
    void printEntries(PrintStream& os, const std::vector<real_t>& m, number_t vb, const SymType sym) const
      {printEntries(os.currentStream(), m, vb, sym);}
    void printEntries(PrintStream& os, const std::vector<complex_t>& m, number_t vb, const SymType sym) const
      {printEntries(os.currentStream(), m, vb, sym);}
    void printEntries(PrintStream& os, const std::vector<Matrix<real_t> >& m, number_t vb, const SymType sym) const
      {printEntries(os.currentStream(), m, vb, sym);}
    void printEntries(PrintStream& os, const std::vector<Matrix<complex_t> >& m, number_t vb, const SymType sym) const
      {printEntries(os.currentStream(), m, vb, sym);}
    void printEntries(PrintStream& os, const std::vector<Vector<real_t> >& m, number_t vb, const SymType sym) const
      {printEntries(os.currentStream(), m, vb, sym);}
    void printEntries(PrintStream& os, const std::vector<Vector<complex_t> >& m, number_t vb, const SymType sym) const
      {printEntries(os.currentStream(), m, vb, sym);}
    void printCooMatrix(PrintStream& os, const std::vector<real_t>& m, SymType s = _noSymmetry) const
      {printCooMatrix(os.currentStream(), m, s);}
    void printCooMatrix(PrintStream& os, const std::vector<complex_t>& m,SymType s = _noSymmetry) const
      {printCooMatrix(os.currentStream(), m, s);}
    void printCooMatrix(PrintStream& os, const std::vector< Matrix<real_t> >& m, SymType s = _noSymmetry) const
      {printCooMatrix(os.currentStream(), m, s);}
    void printCooMatrix(PrintStream& os, const std::vector< Matrix<complex_t> >& m, SymType s = _noSymmetry) const
      {printCooMatrix(os.currentStream(), m, s);}

    // load from file utilities
    //! load from dense matrix in file
    void loadFromFileDense(std::istream& ifs, std::vector<real_t>& mat, SymType sym, bool realAsCmplx)
    {loadSkylineFromFileDense(ifs, mat, rowPointer_, colPointer_, sym, realAsCmplx);}
    //! load from dense matrix in file
    void loadFromFileDense(std::istream& ifs, std::vector<complex_t>& mat, SymType sym, bool realAsCmplx)
    {loadSkylineFromFileDense(ifs, mat, rowPointer_, colPointer_, sym, realAsCmplx);}
    //! load from coordinate matrix in file
    void loadFromFileCoo(std::istream& ifs, std::vector<real_t>& mat, SymType sym, bool realAsCmplx)
    {loadSkylineFromFileCoo(ifs, mat, rowPointer_, colPointer_, sym, realAsCmplx);}
    //! load from coordinate matrix in fire
    void loadFromFileCoo(std::istream& ifs, std::vector<complex_t>& mat, SymType sym, bool realAsCmplx)
    {loadSkylineFromFileCoo(ifs, mat, rowPointer_, colPointer_, sym, realAsCmplx);}


    /*------------------------------------------------------------------------------------------------------------------
                                         basic algebra stuff
      ------------------------------------------------------------------------------------------------------------------*/
    //! Set the Diagonal with a specific value
    template<typename T>
    void setDiagValueDualSkyline(std::vector<T>& v, const T m)
    {
      typename std::vector<T>::iterator it = v.begin() + 1;
      number_t diagSize = diagonalSize() + 1;
      for (number_t k = 1; k < diagSize; k++, it++) {*it = m;}
    }
    //@{
    //! Set value of Diagonal
    void setDiagValue(std::vector<real_t>& m, const real_t k)
    { setDiagValueDualSkyline<>(m, k);}
    void setDiagValue(std::vector<complex_t>& m, const complex_t k)
    { setDiagValueDualSkyline<>(m, k);}
    //@}
    // template Matrix + Matrix addition and specializations (see below)
    template<typename M1, typename M2, typename R>
    void addMatrixMatrix(const std::vector<M1>&, const std::vector<M2>&, std::vector<R>&) const; //!< templated skyline dual Matrix + Matrix

    //@{
    //! Matrix+Matrix
    void addMatrixMatrix(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv) const
    { addMatrixMatrix<>(m, v, rv);}
    void addMatrixMatrix(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv) const
    { addMatrixMatrix<>(m, v, rv); }
    void addMatrixMatrix(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv) const
    { addMatrixMatrix<>(m, v, rv); }
    void addMatrixMatrix(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv) const
    { addMatrixMatrix<>(m, v, rv); }

    template<typename M1, typename M2>
    void addTwoMatrixDualSkyline(std::vector<M1>&, SymType, const std::vector<number_t>&, const std::vector<number_t>&, const std::vector<M2>&, SymType);
    void addTwoMatrix(std::vector<real_t>& m1, SymType st1, const std::vector<number_t>& rowPtr2, const std::vector<number_t>& colPtr2, const std::vector<real_t>& m2, SymType st2)
    { addTwoMatrixDualSkyline<>(m1, st1, rowPtr2, colPtr2, m2, st2); }
    void addTwoMatrix(std::vector<complex_t>& m1, SymType st1, const std::vector<number_t>& rowPtr2, const std::vector<number_t>& colPtr2, const std::vector<real_t>& m2, SymType st2)
    { addTwoMatrixDualSkyline<>(m1, st1, rowPtr2, colPtr2, m2, st2); }
    void addTwoMatrix(std::vector<complex_t>& m1, SymType st1, const std::vector<number_t>& rowPtr2, const std::vector<number_t>& colPtr2, const std::vector<complex_t>& m2, SymType st2)
    { addTwoMatrixDualSkyline<>(m1, st1, rowPtr2, colPtr2, m2, st2); }
    //@}

   /*------------------------------------------------------------------------------------------------------------------
             template Matrix x Vector & Vector x Matrix multiplications and specializations
    ------------------------------------------------------------------------------------------------------------------*/
    template<typename M, typename V, typename R>
    void multMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; //!< templated dual_skyline Matrix x Vector
    template<typename M, typename V, typename R>
    void multVectorMatrix(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; //!< templated Vector x dual_skyline Matrix
    template<typename M, typename V, typename R>
    void multMatrixVector(const std::vector<M>& m, V* vp, R* rp) const;           //!< templated dual_skyline Matrix x Vector (pointer form)
    template<typename M, typename V, typename R>
    void multVectorMatrix(const std::vector<M>& m, V* vp, R* rp) const;           //!< templated Vector x dual_skyline Matrix (pointer form)

    //@{
    //! Matrix (Scalar) * Vector (Scalar)
    void multMatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType sym) const
    { multMatrixVector<>(m, v, rv); }
    void multMatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType sym) const
    { multMatrixVector<>(m, v, rv); }
    void multMatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType sym) const
    { multMatrixVector<>(m, v, rv); }
    void multMatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType sym) const
    { multMatrixVector<>(m, v, rv); }
    //@}

    //@{
    //! Matrix (Matrix) * Vector (Vector)
    void multMatrixVector(const std::vector< Matrix<real_t> >& m, const std::vector< Vector<real_t> >& v, std::vector< Vector<real_t> >& rv, SymType sym) const
    { multMatrixVector<>(m, v, rv); }
    void multMatrixVector(const std::vector< Matrix<real_t> >& m, const std::vector< Vector<complex_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const
    { multMatrixVector<>(m, v, rv); }
    void multMatrixVector(const std::vector< Matrix<complex_t> >& m, const std::vector< Vector<real_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const
    { multMatrixVector<>(m, v, rv); }
    void multMatrixVector(const std::vector< Matrix<complex_t> >& m, const std::vector< Vector<complex_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const
    { multMatrixVector<>(m, v, rv); }
    //@}

    //@{
    //! Vector (Scalar) * Matrix (Scalar)
    void multVectorMatrix(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType sym) const
    { multVectorMatrix<>(m, v, rv); }
    void multVectorMatrix(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType sym) const
    { multVectorMatrix<>(m, v, rv); }
    void multVectorMatrix(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType sym) const
    { multVectorMatrix<>(m, v, rv); }
    void multVectorMatrix(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType sym) const
    { multVectorMatrix<>(m, v, rv); }
    //@}

    //@{ Matrix (Matrix) * Vector (Vector)
    void multVectorMatrix(const std::vector< Matrix<real_t> >& m, const std::vector< Vector<real_t> >& v, std::vector< Vector<real_t> >& rv, SymType sym) const
    { multVectorMatrix<>(m, v, rv); }
    void multVectorMatrix(const std::vector< Matrix<real_t> >& m, const std::vector< Vector<complex_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const
    { multVectorMatrix<>(m, v, rv); }
    void multVectorMatrix(const std::vector< Matrix<complex_t> >& m, const std::vector< Vector<real_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const
    { multVectorMatrix<>(m, v, rv); }
    void multVectorMatrix(const std::vector< Matrix<complex_t> >& m, const std::vector< Vector<complex_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const
    { multVectorMatrix<>(m, v, rv); }
    //@}

    //@{
    //! Matrix * Vector (pointer form)
    void multMatrixVector(const std::vector<real_t>& m, real_t* vp, real_t* rp, SymType sym) const
    { multMatrixVector<>(m, vp, rp); }
    void multMatrixVector(const std::vector<real_t>& m, complex_t* vp, complex_t* rp, SymType sym) const
    { multMatrixVector<>(m, vp, rp); }
    void multMatrixVector(const std::vector<complex_t>& m, real_t* vp, complex_t* rp, SymType sym) const
    { multMatrixVector<>(m, vp, rp); }
    void multMatrixVector(const std::vector<complex_t>& m, complex_t* vp, complex_t* rp, SymType sym) const
    { multMatrixVector<>(m, vp, rp); }
    //@}

     //@{
    //! Vector * Matrix (pointer form)
    void multVectorMatrix(const std::vector<real_t>& m, real_t* vp, real_t* rp, SymType sym) const
    { multVectorMatrix<>(m, vp, rp); }
    void multVectorMatrix(const std::vector<real_t>& m, complex_t* vp, complex_t* rp, SymType sym) const
    { multVectorMatrix<>(m, vp, rp); }
    void multVectorMatrix(const std::vector<complex_t>& m, real_t* vp, complex_t* rp, SymType sym) const
    { multVectorMatrix<>(m, vp, rp); }
    void multVectorMatrix(const std::vector<complex_t>& m, complex_t* vp, complex_t* rp, SymType sym) const
    { multVectorMatrix<>(m, vp, rp); }
    //@}

    /*---------------------------------------------------------------------------------------------------------------------
                                       triangular part matrix * vector
      ------------------------------------------------------------------------------------------------------------------*/
    template<typename M, typename V, typename R>
    void lowerMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
    //@{
    //! lower Matrix (Scalar) * Vector (Scalar)
    virtual void lowerMatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const
    {lowerMatrixVector<>(m,v,rv);}
    virtual void lowerMatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const
    {lowerMatrixVector<>(m,v,rv);}
    virtual void lowerMatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const
    {lowerMatrixVector<>(m,v,rv);}
    virtual void lowerMatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const
    {lowerMatrixVector<>(m,v,rv);}
    //@}
    template<typename M, typename V, typename R>
    void lowerD1MatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
    //@{
    //! lower Matrix diag 1 (Scalar) * Vector (Scalar)
    virtual void lowerD1MatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType s) const
    {lowerD1MatrixVector<>(m,v,rv);}
    virtual void lowerD1MatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType s) const
    {lowerD1MatrixVector<>(m,v,rv);}
    virtual void lowerD1MatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType s) const
    {lowerD1MatrixVector<>(m,v,rv);}
    virtual void lowerD1MatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType s) const
    {lowerD1MatrixVector<>(m,v,rv);}
    //@}
    template<typename M, typename V, typename R>
    void upperMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
    //@{
    //! upper Matrix (Scalar) * Vector (Scalar)
    virtual void upperMatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const
    {upperMatrixVector<>(m,v,rv);}
    virtual void upperMatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const
    {upperMatrixVector<>(m,v,rv);}
    virtual void upperMatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const
    {upperMatrixVector<>(m,v,rv);}
    virtual void upperMatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const
    {upperMatrixVector<>(m,v,rv);}
    //@}
    template<typename M, typename V, typename R>
    void upperD1MatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
    //@{
    //! upper Matrix diag 1 (Scalar) * Vector (Scalar)
   virtual void upperD1MatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const
    {upperD1MatrixVector<>(m,v,rv);}
    virtual void upperD1MatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const
    {upperD1MatrixVector<>(m,v,rv);}
    virtual void upperD1MatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const
    {upperD1MatrixVector<>(m,v,rv);}
    virtual void upperD1MatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const
    {upperD1MatrixVector<>(m,v,rv);}
    //@}
    template<typename M, typename V, typename R>
    void diagonalMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
    //@{
    //! diag Matrix * Vector (Scalar)
    virtual void diagonalMatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const
    {diagonalMatrixVector<>(m,v,rv);}
    virtual void diagonalMatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const
    {diagonalMatrixVector<>(m,v,rv);}
    virtual void diagonalMatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const
    {diagonalMatrixVector<>(m,v,rv);}
    virtual void diagonalMatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const
    {diagonalMatrixVector<>(m,v,rv);}
    //@}

    //----------------------------------------------------------------------------------------------------------------
    //                                     direct solver stuff
    //----------------------------------------------------------------------------------------------------------------
    //! Template LU factorization
    template<typename M>
    void lu(std::vector<M>& m, std::vector<M>& fa) const;

    //! Matrix factorization with LU
    template<typename M>
    void luParallel(std::vector<M>& m, std::vector<M>& fa) const;

    //@{
    //! specializations of template LU
    void lu(std::vector<real_t>& m, std::vector<real_t>& fa, const SymType sym = _noSymmetry) const
    {
       #ifdef XLIFEPP_WITH_OMP
        luParallel<>(m, fa);
       #else
        lu<>(m, fa);
       #endif
    }

    void lu(std::vector<complex_t>& m, std::vector<complex_t>& fa, const SymType sym = _noSymmetry) const
    {
       #ifdef XLIFEPP_WITH_OMP
        luParallel<>(m, fa);
       #else
        lu<>(m, fa);
       #endif
    }
    //@}

    /*!
      Template diagonal and triangular part solvers (after factorization)
      Lower triangular with unit diagonal linear system solver : (I+L) x = b
    */
    template<typename M, typename V, typename X>
    void lowerD1Solver(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
    { lowerD1Solver<>(m, v, x); }
    void lowerD1Solver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { lowerD1Solver<>(m, v, x); }
    void lowerD1Solver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x) const
    { lowerD1Solver<>(m, v, x); }
    void lowerD1Solver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { lowerD1Solver<>(m, v, x); }
    //@}

    //! Diagonal linear system solver : D x = b
    template<typename M, typename V, typename X>
    void diagonalSolver(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
    { diagonalSolver<>(m, v, x); }
    void diagonalSolver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { diagonalSolver<>(m, v, x); }
    void diagonalSolver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x) const
    { diagonalSolver<>(m, v, x); }
    void diagonalSolver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const
    { diagonalSolver<>(m, v, x); }
    //@}

    //! Upper triangular with unit diagonal linear system solver : (I+U) x = b
    template<typename M, typename V, typename X>
    void upperD1Solver(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym) 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
    { upperD1Solver<>(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
    { upperD1Solver<>(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
    { upperD1Solver<>(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
    { upperD1Solver<>(m, v, x, sym); }
    //@}

    //! upper triangular linear system solver : (D+U) x = b
    template<typename M, typename V, typename X>
    void upperSolver(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) 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) const
    { upperSolver<>(m, v, x); }
    void upperSolver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType) const
    { upperSolver<>(m, v, x); }
    void upperSolver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x, const SymType sym) const
    { upperSolver<>(m, v, x); }
    void upperSolver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym) const
    { upperSolver<>(m, v, x); }
    //@}

    /* ----------------------------------------------------------------------------------------------------------------------------
                                                 UMFpack stuff
       ----------------------------------------------------------------------------------------------------------------------------*/
    //! conversion to umfpack format
    template<typename M, typename OrdinalType>
    void toUmfPack(const std::vector<M>& values, std::vector<OrdinalType>& colPointer, std::vector<OrdinalType>& rowIndex, std::vector<M>& mat) const;
    //@{
    //! specializations of umfpack conversion
    void toUmfPack(const std::vector<real_t>& values, std::vector<int_t>& colPointer, std::vector<int_t>& rowIndex, std::vector<real_t>& mat, const SymType sym) const
    { toUmfPack<>(values,colPointer,rowIndex,mat);}
    void toUmfPack(const std::vector<complex_t>& values, std::vector<int_t>& colPointer, std::vector<int_t>& rowIndex, std::vector<complex_t>& mat, const SymType sym) const
    { toUmfPack<>(values,colPointer,rowIndex,mat);}
    //@}

}; //end of class definition


// constructor from the list of col index of non zero for each row
//   cols[k] of type L stores the indices of column of non zero element on row k (index starts at 1!)
//   only requirement on L are the definition of forward const iterator named ::const_iterator
//                             and the size
//   works for L = vector<number_t>, list<number_t>, set<number_t>
template<class L>
DualSkylineStorage::DualSkylineStorage(number_t nr, number_t nc, const std::vector<L>& cols, string_t id)
  : SkylineStorage(nr, nc, _dual, id)
{
  trace_p->push("DualSkylineStorage constructor");

  // construct rowPointer_ for strict lower part
  typename std::vector<L>::const_iterator itvn;
  rowPointer_.resize(nbRows_ + 1);
  colPointer_.resize(nbCols_ + 1);
  std::vector<number_t>::iterator it = rowPointer_.begin();
  number_t r = 2, l = 0;
  *it = 0; it++;
  for(itvn = cols.begin()+1; itvn != cols.end(); itvn++, r++, it++) {
      *it = *(it - 1) + l;
      l=0;
      if(!itvn->empty()) { // length
            number_t s = *std::min_element(itvn->begin(), itvn->end());
            if(s<r) l = r-s;
      }
    }
  *it = *(it - 1) + l;

  //construct colPointer_ for strict upper part
  number_t c = 1;
  for(it = colPointer_.begin(); it != colPointer_.end(); it++,c++) { *it = c; }
  it=colPointer_.begin();
  typename L::const_iterator its;
  r=1;
  for(itvn = cols.begin(); itvn != cols.end(); itvn++, r++)
    {
      for(its = itvn->begin(); its != itvn->end(); its++)
        if(*its > r) //strict upper part
          {
            c=*its-1;
            *(it+c) = std::min(*(it+c), r); // update min row index on col c
          }
    }

  *(colPointer_.begin()) = 0;
  c = 2; l=0;
  for(it = colPointer_.begin() + 1; it != colPointer_.end(); it++, c++)
    {
      r=*it;
      *it = *(it - 1) + l;
      l=c-r;
    }

  trace_p->pop();
}

} // end of namespace xlifepp

#endif // DUAL_SKYLINE_STORAGE_HPP