xref: /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/largeMatrix/LargeMatrix.hpp

/*
XLiFE++ issue 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 LargeMatrix.hpp
  \author E. Lunéville, M.-H. Nguyen
  \since 24 jun 2012
  \date 6 Juin 2014

  \brief Definition of the xlifepp::LargeMatrix class

  Class to deal with large matrices. This class mainly handles a vector of values
  and a pointer to a storage (see storage classes). Algorithms are defined in storage classes
  The class is templated by the type of the values real/complex scalar/vector/matrix

  IMPORTANT NOTE : by convenience, the first adress of values_ vector (say values_[0]) always exists and contains the 0 value of type T
  the actual matrix values begins at adress 1
*/

#ifndef LARGE_MATRIX_HPP
#define LARGE_MATRIX_HPP

#include "config.h"
#ifdef XLIFEPP_WITH_UMFPACK
#include "umfpackSupport.h"
#endif
#include "MatrixStorage.hpp"
#include "denseStorage/DenseStorage.hpp"
#include "denseStorage/ColDenseStorage.hpp"
#include "denseStorage/RowDenseStorage.hpp"
#include "denseStorage/DualDenseStorage.hpp"
#include "denseStorage/SymDenseStorage.hpp"
#include "csStorage/CsStorage.hpp"
#include "csStorage/ColCsStorage.hpp"
#include "csStorage/RowCsStorage.hpp"
#include "csStorage/DualCsStorage.hpp"
#include "csStorage/SymCsStorage.hpp"
#include "skylineStorage/DualSkylineStorage.hpp"
#include "skylineStorage/SymSkylineStorage.hpp"
#include "eigenSolvers.h"
#include "eigenSparseInterface/LargeMatrixAdapter.hpp"
#include "eigenSparseInterface/LargeMatrixAdapterInverse.hpp"
#include "eigenSparseInterface/LargeMatrixAdapterInverseGeneralized.hpp"
#include "eigenSparseInterface/MultiVectorAdapter.hpp"
#include "../arpackppSupport/arpackSolve.hpp"
#include "utils.h"

namespace xlifepp
{

/*!
  \struct LargeMatrixTrait
  A small struct to deal with type of scalar (real, complex) of LargeMatrix
*/
template<typename K>
struct LargeMatrixTrait
{};

#ifndef DOXYGEN_SHOULD_SKIP_THIS

template<>
struct LargeMatrixTrait<real_t>
{
  typedef real_t ScalarType;
};

template<>
struct LargeMatrixTrait<complex_t>
{
  typedef complex_t ScalarType;
};

template<>
struct LargeMatrixTrait<Matrix<real_t> >
{
  typedef real_t ScalarType;
};

template<>
struct LargeMatrixTrait<Matrix<complex_t> >
{
  typedef complex_t ScalarType;
};

#endif

//-----------------------------------------------------------------------------------------
/*!
   \class LargeMatrix
   templated class to deal with large matrices,
   it handles a vector of values (non zero coefficients) and a pointer to a MatrixStorage
 */
//-----------------------------------------------------------------------------------------
template <typename T>
class LargeMatrix
{
public:
  ValueType valueType_;     //!< type of values (real, complex)
  StrucType strucType;      //!< struc of values (scalar, vector, matrix)
  number_t nbRows;          //!< number of rows counted in T
  number_t nbCols;          //!< number of columns counted in T
  SymType sym;              //!< type of symmetry
  dimen_t nbRowsSub;        //!< number of rows of submatrices (1 if scalar values)
  dimen_t nbColsSub;        //!< number of columns cof submatrices (1 if scalar values)
  string_t name;            //!< optionnal name, useful for documentation

  //for factorized form
  FactorizationType factorization_;     //!< one of _noFactorization, _lu, _ldlt, _ldlstar, _llt, _llstar, _qr, _umfpack
  std::vector<number_t> rowPermutation_;//!< row permutation given by some factorizations (LU in row dense storage for instance)
  std::vector<number_t> colPermutation_;//!< col permutation given by some factorizations (LU in col dense storage for instance)
#ifdef XLIFEPP_WITH_UMFPACK
  UmfPackLU<LargeMatrix<T> >* umfLU_p;  //!< pointer to store LU factors built by UmfPack
#endif

  typedef typename LargeMatrixTrait<T>::ScalarType ScalarType; //!< useful typedef for LargeMatrix class

protected:
  std::vector<T> values_;   //!< list of values ordered along storage choice
  MatrixStorage* storage_p; //!< pointer to the matrix storage
  void allocateStorage(StorageType sto, AccessType at, const T& val = T());  //!< allocate storage if not exists
  void setType(const T&);                                                    //!< initialize type of matrix elements

public:
  // constructors, destructor
  LargeMatrix();                                                             //!< default constructor
  LargeMatrix(ValueType vt,StrucType st, number_t nbr, number_t nbc, SymType sy,
              dimen_t nbrs, dimen_t nbcs, const T& val, const string_t& na, MatrixStorage* sp=0);  //!< explicit constructor
  //LargeMatrix(const LargeMatrix<T>&);                                      //!< copy constructor (shadow copy)
  LargeMatrix(const LargeMatrix<T>&, bool storageCopy = false);              //!< copy constructor (shadow copy), with option to copy or not the storage
  template<typename K> LargeMatrix(const LargeMatrix<K>&, bool storageCopy = false); //!< generalized copy constructor (shallow copy), with option to copy or not the storage

  LargeMatrix(MatrixStorage*, const T& val = T(0), SymType sy = _noSymmetry);//!< construct matrix from storage and constant value
  LargeMatrix(MatrixStorage* ms, dimPair dims, SymType sy = _noSymmetry);    //!< construct matrix from storage and symmetry type (matrix values)
  LargeMatrix(MatrixStorage*, SymType);                                      //!< construct matrix from storage and symmetry type (scalar values assumed)
  ~LargeMatrix();                                                            //!< destructor
  //! construct matrix from dimensions, storage type and a constant value
  LargeMatrix(number_t nr, number_t nc, StorageType sto = _dense, AccessType at = _row, const T& val = T());
  //! construct "symmetric" matrix from dimensions, storage type, symmetry type a constant value
  LargeMatrix(number_t nr, number_t nc, StorageType sto = _dense, SymType = _symmetric, const T& val = T());
  //! construct matrix from file
  LargeMatrix(const string_t& fn, StorageType fst, number_t nr, number_t nc, StorageType sto = _dense, AccessType at = _row, number_t nsr = 1, number_t nsc = 1);
  //! construct "symmetric" matrix from file
  LargeMatrix(const string_t& fn, StorageType fst, number_t nr, StorageType sto = _dense, SymType sy = _symmetric, number_t nsr = 1);
  //! construct Special matrix (diagonal)
  LargeMatrix(SpecialMatrix sp, StorageType st, AccessType at, number_t nbr, number_t nbc, const T v);

  void init(MatrixStorage* ms, const T& val, SymType sy); //!< initialize matrix from storage
  void clearall();         //!< clear values and storage pointer (not the attributes)
  void clear()             //! deallocate memory used by matrix values, storage is not deallocated!
  {
    if(Trace::traceMemory)
      {
        thePrintStream<<"LargeMatrix::clear de-allocates a large matrix : "<<&values_<<", "
                      <<values_.size()<<" non zeros coefficients "<<dimValues();
        if(storage_p!=0) thePrintStream<<", storage "<<storage_p->name();
        thePrintStream<<eol<<std::flush;
      }
    std::vector<T> empty;
    values_.swap(empty);
  }

  LargeMatrix<T>& operator=(const LargeMatrix<T>&); //!< assign operator
  template<typename K> LargeMatrix<T>& operator=(const LargeMatrix<K>&); //!< generalized assign operator

  // various utilities
  ValueType valueType() const                //! return type of values (real, complex)
  {
    return valueType_;
  }

  std::vector<T>& values()                   //! access to values
  {
    return values_;
  }

  const std::vector<T>& values() const       //! access to values (const)
  {
    return values_;
  }

  number_t numberOfRows() const              //! return number of rows counted in T
  {
    return nbRows;
  }

  number_t numberOfCols() const              //! return number of columns counted in T
  {
    return nbCols;
  }

  SymType symmetry() const                   //! return kind of symmetry
  {
    return sym;
  }

  dimPair dimValues() const                  //! return dimensions of values, (1,1) when scalar
  {
    return dimPair(nbRowsSub, nbColsSub);
  }

  StorageType storageType() const            //! return the storage type
  {
    return storage_p->storageType();
  }

  AccessType accessType() const              //! return the access type
  {
    return storage_p->accessType();
  }

  MatrixStorage* storagep() const            //! return the storage pointer (const)
  {
    return storage_p;
  }

  MatrixStorage*& storagep()                 //! return the storage pointer
  {
    return storage_p;
  }

  number_t nbNonZero() const                 //! number of non zero coefficients stored (given by storage)
  {
    return values_.size() - 1;
  }

  void resetZero()                           //! reset first value to 0
  {
    values_[0] *= 0.;
  }

  //access to values or addresses of values
  T operator()(number_t i, number_t j) const //! access to entry (i,j) (do not use in heavy computation)
  {
    return get(i,j);
  }

  T get(number_t i, number_t j) const //! access to entry (i,j) (do not use in heavy computation)
  {
    if(sym!=_noSymmetry && storage_p->accessType()==_sym)
      switch(sym)
        {
        case _skewSymmetric :
          return -values_[storage_p->pos(i, j, sym)];
          break;
        case _selfAdjoint   :
          return conj(values_[storage_p->pos(i, j, sym)]);
          break;
        case _skewAdjoint   :
          return -conj(values_[storage_p->pos(i, j, sym)]);
          break;
        default :
          return values_[storage_p->pos(i, j, sym)];
        }
    return values_[storage_p->pos(i, j, sym)];
  }

  T& operator()(number_t i, number_t j, bool errorOn=false) //! access to entry (i,j) with check (do not use in heavy computation)
  {
    number_t posij = storage_p->pos(i, j, sym);
    if(errorOn)
      if(posij == 0 || (sym!=_noSymmetry && storage_p->accessType()==_sym && j>i))
        {
          error("largematrix_indicesout", name, i, j);
        }
    return values_[posij];
  }

  number_t pos(number_t i, number_t j) const
  {
    return storage_p->pos(i, j, sym);
  }

  std::vector<std::pair<number_t, number_t> > getCol(number_t c, number_t r1=1, number_t r2=0) const //! get (row indices, adresses) of col c in set [r1,r2] (r2=0 means nbOfRows)
  {
    return storage_p->getCol(sym,c,r1,r2);
  }

  std::vector<std::pair<number_t, number_t> > getRow(number_t r, number_t c1=1, number_t c2=0) const //! /! get (col indices, adresses) of row r in set [c1,c2] (c2=0 means nbOfCols)
  {
    return storage_p->getRow(sym,r,c1,c2);
  }

  std::vector<T> col(number_t c) const;       //!< get col c (>=1) as vector
  std::vector<T> row(number_t r) const;       //!< get row r (>=1) as vector

  T operator()(number_t adr) const              //! access to entry at adress adr
  {
    return values_[adr];
  }

  T& operator()(number_t adr)                   //! access to entry at adress adr
  {
    return values_[adr];
  }

  typename std::vector<T>::iterator at(number_t i, number_t j)
  {
    number_t p=pos(i,j);
    if(p>0) return values_.begin()+ p;
    else return values_.end();
  }

  typename std::vector<T>::const_iterator at(number_t i, number_t j) const
  {
    number_t p=pos(i,j);
    if(p>0) return values_.begin()+ p;
    else return values_.end();
  }

  // give the positions in storage of a submatrix given by its row and column index
  // the matrix of positions (adrs) are returned in a vector (row matrix storage)
  // as xlifepp::Matrix inherits from std::vector, it can be used as input argument
  // if errorOn is true, indices out of storage raises an error
  void positions(const std::vector<number_t>& rows, const std::vector<number_t>& cols,
                 std::vector<number_t>& adrs, bool errorOn=false) const //! positions in storage of a submatrix
  {
    storage_p->positions(rows, cols, adrs, errorOn, sym);
  }

  //removing operation
  void deleteRows(number_t, number_t);             //!< delete rows r1,..., r2 (may be expansive for some storage)
  void deleteCols(number_t, number_t);             //!< delete cols c1,..., c2 (may be expansive for some storage)
  void setColToZero(number_t c1=0, number_t c2=0); //!< set to 0 cols c1,..., c2 (index start from 1), no col deletion
  void setRowToZero(number_t r1=0, number_t r2=0); //!< set to 0 rows r1,..., r2 (index start from 1), no row deletion

  //storage conversion
  void toSkyline();                         //!< convert current matrix storage to skyline storage
  void toStorage(MatrixStorage*);           //!< convert current matrix storage to an other matrix storage
  template<typename K>
  LargeMatrix<K>* toScalar(K);              //!< create new scalar matrix entry from non scalar matrix entry
  template <typename K>
  void toScalarEntries(const std::vector<Matrix<K> >&, std::vector<K>&, const MatrixStorage&); //!< copy matrix values to scalar values
  void toUnsymmetric(AccessType at=_sym);      //!< convert current matrix to unsymmetric representation if it has a symmetric one
  bool isDiagonal() const;                     //!< true if matrix is in fact diagonal
  bool isId(const double & tol=0.) const;      //!< true if matrix is in fact id with a tolerance
  void roundToZero(real_t aszero=10*theEpsilon); //!< round to 0 all coefficients smaller (in norm) than a aszero, storage is not modified
  LargeMatrix<T>& toConj();                      //!< conjugate matrix

  void extract2UmfPack(std::vector<int_t>& colPointer, std::vector<int_t>& rowIndex, std::vector<T>& matA) const;   //!< extract and convert largeMatrix (and its storage) to Umfpack
  bool getCsColUmfPack(std::vector<int_t>& colPointer, std::vector<int_t>& rowIndex, const T*& matA) const;         //!< get largeMatrix and its storage to Umfpack (only for cs col)

  LargeMatrix<T>* extract(const std::vector<number_t>& rowIndex,
                          const std::vector<number_t>& colIndex) const;              //!< extract LargeMatrix from current LargeMatrix (return a pointer to LargeMatrix)
  LargeMatrix<T>* extract(number_t r1, number_t r2, number_t c1, number_t c2) const; //!< extract LargeMatrix from current LargeMatrix (return a pointer to LargeMatrix)

  // algebraic utilities
  template<typename iterator>
  LargeMatrix<T>& add(iterator&,iterator&);                            //!< add values given by iterator from beginning
  LargeMatrix<T>& add(const std::vector<T>&, const std::vector<number_t>&,
                      const std::vector<number_t>&);                     //!< add a matrix at ranks (rows, cols)
  LargeMatrix<T>& add(const T& c, const std::vector<number_t>& rows,
                      const std::vector<number_t>& cols);                //!< add a constant to submatrix given by its ranks (rows, cols)
  template <typename K, typename C>
  LargeMatrix<T>& add(const LargeMatrix<K>&,
                      const std::vector<number_t>&,
                      const std::vector<number_t>&, C);                  //!< add scaled LargeMatrix at ranks (rows, cols) of current LargeMatrix
  void addColToCol(number_t c1, number_t c2, complex_t a, bool updateStorage=false); //!< combine cols of matrix : c2 = c2 + a *c1 (indices start from 1) with update storage option
  void addRowToRow(number_t r1, number_t r2, complex_t a, bool updateStorage=false); //!< combine rows of matrix : r2 = r2 + a *r1 (indices start from 1) with update storage option
  void copyVal(const LargeMatrix<T>&, const std::vector<number_t>&,
               const std::vector<number_t>&);                                    //!< copy in current matrix values of mat located at mat_adrs at position cur_adrs
  template <typename K>
  LargeMatrix<T>& assign(const LargeMatrix<K>&, const std::vector<number_t>&, const std::vector<number_t>&); //!< assign LargeMatrix at ranks (rows, cols) to current LargeMatrix

  template<typename K> LargeMatrix<T>& operator*=(const K&);           //!< LargeMatrix<T> *= v
  template<typename K> LargeMatrix<T>& operator/=(const K&);           //!< LargeMatrix<T> /= v
  LargeMatrix<T>& operator+=(const LargeMatrix<T>&);                   //!< LargeMatrix<T> += LargeMatrix<T>
  LargeMatrix<T>& operator-=(const LargeMatrix<T>&);                   //!< LargeMatrix<T> -= LargeMatrix<T>

  real_t norm2() const;                                                //!< Frobenius norm of a LargeMatrix
  real_t norminfty() const;                                            //!< infinite norm of a LargeMatrix
  real_t partialNormOfCol(number_t c, number_t r1, number_t r2) const; //!< partial norm of column c, only components in range [r1,r2]
  real_t squaredNorm() const                                           //! square of Frobenius norm of a LargeMatrix
  {
    real_t n=norm2();
    return n*n;
  }

  void ldltFactorize();                        //!< Factorization LDLt
  void ldlstarFactorize();                     //!< Factorization LDL*
  void luFactorize(bool withPermutation=true); //!< Factorization LU
  void iluFactorize();                         //!< Incomplete Factorization LU
  void ildltFactorize();                       //!< Incomplete Factorization LDLt
  void illtFactorize();                        //!< Incomplete Factorization LLt
  void ildlstarFactorize();                    //!< Incomplete Factorization LDL*
  void illstarFactorize();                     //!< Incomplete Factorization LL*
  void umfpackFactorize();                     //!< Factorization using umfpack

  void multMatrixCol(T* pM, T* pR, number_t p) const;     //!< product LargeMatrix * dense col matrix passed as pointer
  void multLeftMatrixCol(T* pM, T* pR, number_t p) const; //!< product dense col matrix * LargeMatrix passed as pointer
  void multMatrixRow(T* pM, T* pR, number_t p) const;     //!< product LargeMatrix * dense row matrix passed as pointer
  void multLeftMatrixRow(T* pM, T* pR, number_t p) const; //!< product dense row matrix * LargeMatrix passed as pointer

  // extern templated product matrix*vector and "specializations" real_t/complex_t and Scalar/Vector/Matrix
  template <typename S, typename V, typename R>
  friend void multMatrixVector(const LargeMatrix<S>& mat, V*, R*);
  template <typename S>
  friend void multMatrixVector(const LargeMatrix<S>& mat, const std::vector<S>& vec, std::vector<S>& res);
  template <typename S>
  friend void multMatrixVector(const LargeMatrix<Matrix<S> >& mat, const std::vector<Vector<S> >& vec, std::vector<Vector<S> >& res);
  friend void multMatrixVector(const LargeMatrix<complex_t>& mat, const std::vector<real_t>& vec, std::vector<complex_t>& res);
  friend void multMatrixVector(const LargeMatrix<real_t>& mat, const std::vector<complex_t>& vec, std::vector<complex_t>& res);
  template <typename S>
  friend void multMatrixVector(const LargeMatrix<Matrix<S> >& mat, const std::vector<Vector<S> >& vec, std::vector<Vector<S> >& res);
  friend void multMatrixVector(const LargeMatrix<Matrix<complex_t> >& mat, const std::vector<Vector<real_t> >& vec, std::vector<Vector<complex_t> >& res);
  friend void multMatrixVector(const LargeMatrix<Matrix<real_t> >& mat, const std::vector<Vector<complex_t> >& vec, std::vector<Vector<complex_t> >& res);

  // extern templated product vector*matrix and "specializations" real_t/complex_t and Scalar/Vector/Matrix
  template <typename S, typename V, typename R>
  friend void multVectorMatrix(V*, const LargeMatrix<S>& mat, R*);
  template <typename S>
  friend void multVectorMatrix(const std::vector<S>& vec, const LargeMatrix<S>& mat, std::vector<S>& res);
  friend void multVectorMatrix(const std::vector<real_t>& vec, const LargeMatrix<complex_t>& mat, std::vector<complex_t>& res);
  friend void multVectorMatrix(const std::vector<complex_t>& vec, const LargeMatrix<real_t>& mat, std::vector<complex_t>& res);
  template <typename S>
  friend void multVectorMatrix(const std::vector<Vector<S> >& vec, const LargeMatrix<Matrix<S> >& mat, std::vector<Vector<S> >& res);
  friend void multVectorMatrix(const std::vector<Vector<real_t> >& vec, const LargeMatrix<Matrix<complex_t> >& mat, std::vector<Vector<complex_t> >& res);
  friend void multVectorMatrix(const std::vector<Vector<complex_t> >& vec, const LargeMatrix<Matrix<real_t> >& mat, std::vector<Vector<complex_t> >& res);
  template <typename S, typename V, typename R>
  friend void multVectorMatrix(const LargeMatrix<S>& mat, V*, R*);
  template <typename S>
  friend void multVectorMatrix(const LargeMatrix<S>& mat, const std::vector<S>& vec, std::vector<S>& res);
  friend void multVectorMatrix(const LargeMatrix<complex_t>& mat, const std::vector<real_t>& vec, std::vector<complex_t>& res);
  friend void multVectorMatrix(const LargeMatrix<real_t>& mat, const std::vector<complex_t>& vec, std::vector<complex_t>& res);
  template <typename S>
  friend void multVectorMatrix(const LargeMatrix<Matrix<S> >& mat, const std::vector<Vector<S> >& vec, std::vector<Vector<S> >& res);
  friend void multVectorMatrix(const LargeMatrix<Matrix<complex_t> >& mat, const std::vector<Vector<real_t> >& vec, std::vector<Vector<complex_t> >& res);
  friend void multVectorMatrix(const LargeMatrix<Matrix<real_t> >& mat, const std::vector<Vector<complex_t> >& vec, std::vector<Vector<complex_t> >& res);

  // LargeMatrix MUST be FACTORIZED before being used
  template <typename S1, typename S2, typename S3>
  friend void multFactMatrixVector(const LargeMatrix<S1>& mat, const std::vector<S2>& vec, std::vector<S3>& res);
  template <typename S1, typename S2, typename S3>
  friend void multVectorFactMatrix(const LargeMatrix<S1>& mat, const std::vector<S2>& vec, std::vector<S3>& res);
  template<typename S1, typename S2>
  friend void multInverMatrixVector(const LargeMatrix<S1>& mat, const std::vector<S2>& vec,
                                    std::vector<typename Conditional<NumTraits<S1>::IsComplex, S1, S2>::type>& res,
                                    FactorizationType fac);

  // extern template product matrix*matrix and "specializations" real_t/complex_t and Scalar/Matrix
  template <typename SA, typename SB, typename SR>
  friend void multMatrixMatrix(const LargeMatrix<SA>&, const LargeMatrix<SB>&, LargeMatrix<SR>&);
  template <typename SA, typename SB, typename SR>
  friend void multMatrixMatrix(const LargeMatrix<SA>&, const std::vector<SB>&, LargeMatrix<SR>&);
  template <typename SA, typename SB, typename SR>
  friend void multMatrixMatrix(const std::vector<SA>&, const LargeMatrix<SB>&, LargeMatrix<SR>&);

  //! templated factorization LDLt solver
  template<typename S1, typename S2>
  void ldltSolve(std::vector<S1>& vec, std::vector<S2>& res) const;
  //! templated factorization LDL* solver
  template<typename S>
  void ldlstarSolve(std::vector<S>& vec, std::vector<T>& res) const;
  //! templated factorization LLt solver
  template<typename S1, typename S2>
  void lltSolve(std::vector<S1>& vec, std::vector<S2>& res) const;
  //! templated factorization LU solver
  template<typename S1, typename S2>
  void luSolve(std::vector<S1>& vec, std::vector<S2>& res) const;
  //! templated factorization LU solver (UmfPack)
  template<typename S1, typename S2>
  void umfluSolve(std::vector<S1>& vec, std::vector<S2>& res) const;
  //! templated factorization SOR lower part solver
  template<typename S1, typename S2>
  void sorLowerMatrixVector(const std::vector<S1>& vec, std::vector<S2>& res, const real_t w) const;
  //! templated factorization SOR diagonal solver
  template<typename S1, typename S2>
  void sorDiagonalMatrixVector(const std::vector<S1>& vec, std::vector<S2>& res, const real_t w) const;
  //! templated factorization SOR upper part solver
  template<typename S1, typename S2>
  void sorUpperMatrixVector(const std::vector<S1>& vec, std::vector<S2>& res, const real_t w) const;
  //! templated factorization SOR lower part solver
  template<typename S1, typename S2>
  void sorLowerSolve(const std::vector<S1>& vec, std::vector<S2>& res, const real_t w) const;
  //! templated factorization SOR diagonal solver
  template<typename S1, typename S2>
  void sorDiagonalSolve(const std::vector<S1>& vec, std::vector<S2>& res, const real_t w) const;
  //! templated factorization SOR upper part solver
  template<typename S1, typename S2>
  void sorUpperSolve(const std::vector<S1>& vec, std::vector<S2>& res, const real_t w) const;

  // extern templated eigen Davidson solver
  template<typename K>
  friend number_t eigenDavidsonSolve(const LargeMatrix<K>* pA, const LargeMatrix<K>* pB, std::vector<std::pair<complex_t,Vector<complex_t> > >& res,
                                     number_t nev, real_t tol, string_t which, bool isInverted, FactorizationType fac, bool isShift);
  // extern templated eigen Krylov-Schur solver
  template<typename K>
  friend number_t eigenKrylovSchurSolve(const LargeMatrix<K>* pA, const LargeMatrix<K>* pB, std::vector<std::pair<complex_t,Vector<complex_t> > >& res,
                                        number_t nev, real_t tol, string_t which, bool isInverted, FactorizationType fac, bool isShift);

  //@{
  //! templated UMFPACK solver
  template<typename S>
  real_t umfpackSolve(const std::vector<S>& vec,
                      std::vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, S>::type>& res,
                      bool soa=true);    //!< solve linear system using umfpack
  template<typename S>
  real_t umfpackSolve(const std::vector<std::vector<S>*>& bs,
                      std::vector<std::vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType,S>::type>* >& xs,
                      bool soa=true);   //!< solve linear system using umfpack, multiple rhs
  template<typename S>
  real_t umfpackSolve(const std::vector<int_t>& colPointer, const std::vector<int_t>& rowIndex,
                      const std::vector<T>& values, const std::vector<S>& vec,
                      std::vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, S>::type>& res,
                      bool soa=true); //!< solve linear system using umfpack
  //@}

  //@{
  //! lower and upper solvers on vector : LD1^{-1}*b, U^{-1}*b, b*LD1^{-1}, b*U^{-1}, b may be modified
  std::vector<T>& lowerD1Solve(std::vector<T>& b, std::vector<T>& x) const;
  std::vector<T>& upperSolve(std::vector<T>& b, std::vector<T>& x) const;
  std::vector<T>& lowerD1LeftSolve(std::vector<T>& b, std::vector<T>& x) const;
  std::vector<T>& upperLeftSolve(std::vector<T>& b, std::vector<T>& x) const;
  //@}

  //@{
  //! lower and upper solvers on LargeMatrix : LD1^{-1}*M, U^{-1}*M, M*LD1^{-1}, M*U^{-1}, result is stored as col dense
  LargeMatrix<T>& lowerD1Solve(const LargeMatrix<T>& M, const LargeMatrix<T>& X) const;
  LargeMatrix<T>& upperSolve(const LargeMatrix<T>& M, const LargeMatrix<T>& X) const;
  LargeMatrix<T>& lowerD1LeftSolve(const LargeMatrix<T>& M, const LargeMatrix<T>& X) const;
  LargeMatrix<T>& upperLeftSolve(const LargeMatrix<T>& M, const LargeMatrix<T>& X) const;
  //@}

  // external template addition two matrices A+B=C, all of them have the same type of storage
  template<typename S>
  friend void addMatrixMatrix(const LargeMatrix<S>& matA, const LargeMatrix<S>& matB, LargeMatrix<S>& matC);
#ifndef DOXYGEN_SHOULD_SKIP_THIS
  friend void addMatrixMatrix(const LargeMatrix<complex_t>& matA, const LargeMatrix<real_t>& matB, LargeMatrix<complex_t>& matC);
  friend void addMatrixMatrix(const LargeMatrix<real_t>& matA, const LargeMatrix<complex_t>& matB, LargeMatrix<complex_t>& matC);
#endif

  // Scale a largeMatrix with a scalar
  template<typename S>
  friend LargeMatrix<S> multMatrixScalar(const LargeMatrix<S>&, const S);
#ifndef DOXYGEN_SHOULD_SKIP_THIS
  friend LargeMatrix<complex_t> multMatrixScalar(const LargeMatrix<complex_t>&, const real_t);
  friend LargeMatrix<complex_t> multMatrixScalar(const LargeMatrix<real_t>&, const complex_t);
#endif

  // Construct a diagonalMatrix from a Matrix
  template<typename S>
  friend LargeMatrix<S> diagonalMatrix(const LargeMatrix<S>&, const S);

  // Construct a identity matrix from a Matrix
  template<typename S>
  friend LargeMatrix<S> identityMatrix(const LargeMatrix<S>&);

  // print utilities
  void print(std::ostream&) const;                                       //!< print on output stream
  void print(PrintStream& os) const
  { print(os.currentStream()); }
  void viewStorage(std::ostream&) const;                                 //!< visualize storage on output stream
  void viewStorage(PrintStream& os) const
  { viewStorage(os.currentStream()); }
  template <typename S>
  friend std::ostream& operator<<(std::ostream&, const LargeMatrix<S>&);

  // save to file utilities
  void saveToFile(const string_t&, StorageType, bool encode = false) const; //!< save matrix to a file along different format
  void saveToFile(const char*, StorageType, bool encode = false) const;   //!< save matrix to a file along different format
  void loadFromFile(const string_t&, StorageType);                          //!< load matrix from file
  string_t encodeFileName(const string_t&, StorageType) const;                //!< encode file name with additional matrix informations
};

// utilities
//std::pair<dimen_t, dimen_t> dimsOf(const real_t&);
//std::pair<dimen_t, dimen_t> dimsOf(const complex_t&);
//std::pair<dimen_t, dimen_t> dimsOf(const Matrix<real_t>&);
//std::pair<dimen_t, dimen_t> dimsOf(const Matrix<complex_t>&);

// external template addition of two matrix A+B=C, they need not share a same storage
// This function return a non-null pointer whose referenced memory should be MANAGEd after being called.
template<typename S>
LargeMatrix<S>* addMatrixMatrixSkyline(const LargeMatrix<S>& matA, const LargeMatrix<S>& matB);

/*---------------------------------------------------------------------------------------------------
  LargeMatrix<T> constructors, destructor
---------------------------------------------------------------------------------------------------*/
//default constructor (no storage)
template <typename T>
LargeMatrix<T>::LargeMatrix()
  : valueType_(_real), strucType(_scalar), nbRows(0), nbCols(0), sym(_noSymmetry), nbRowsSub(1), nbColsSub(1)
{
  setType(T());
  storage_p = 0;
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  factorization_=_noFactorization;
  values_.resize(1, T());
}
//explicit constructor
template <typename T>
LargeMatrix<T>::LargeMatrix(ValueType vt,StrucType st, number_t nbr, number_t nbc, SymType sy,
                            dimen_t nbrs, dimen_t nbcs, const T& val, const string_t& na,  MatrixStorage* sp)
  : valueType_(vt),strucType(st),nbRows(nbr),nbCols(nbc),sym(sy),nbRowsSub(nbrs),nbColsSub(nbcs), name(na), storage_p(sp)
{
  if(sp!=0) init(sp,val,sy);
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  factorization_=_noFactorization;
}

template <typename T>
LargeMatrix<T>::LargeMatrix(const LargeMatrix<T>& mat, bool storageCopy)
  : valueType_(mat.valueType_), strucType(mat.strucType), nbRows(mat.nbRows), nbCols(mat.nbCols), sym(mat.sym),
    nbRowsSub(mat.nbRowsSub), nbColsSub(mat.nbColsSub), name("shallow copy of" + mat.name)
{
  setType(T());
  nbRowsSub=mat.nbRowsSub;  //redo because setType reset to 1
  nbColsSub=mat.nbColsSub;  //redo because setType reset to 1

  if (Trace::traceMemory)
    {
      thePrintStream << "LargeMatrix::copy_constructor allocates a new large matrix : " << &values_ << ", "
                     << mat.values_.size() << " non zeros coefficients " << dimValues();
      if(mat.storage_p != 0) thePrintStream << ", storage " << mat.storage_p->name();
      thePrintStream << eol << std::flush;
    }

  values_ = mat.values_;
  storage_p = mat.storage_p;
  if (storage_p == 0) return;
  if (storageCopy) storage_p=mat.storage_p->clone();               // copy storage
  storage_p->objectPlus();
  factorization_=mat.factorization_;
  #ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
  if (mat.umfLU_p != 0)
    {
        warning("free_warning","LargeMatrix copy, does not copy the umfpackLU object, redo factorisation if required");
        factorization_=_noFactorization;
    }
  #endif
}

template <typename T>
LargeMatrix<T>& LargeMatrix<T>::operator=(const LargeMatrix<T>& mat)
{
  if (this != &mat)
    {
      clearall();
      valueType_=mat.valueType_;
      strucType=mat.strucType;
      nbRows=mat.nbRows;
      nbCols=mat.nbCols;
      sym=mat.sym;
      nbRowsSub=mat.nbRowsSub;
      nbColsSub=mat.nbColsSub;
      name="shallow copy of" + mat.name;
      setType(T());
      nbRowsSub=mat.nbRowsSub;  //redo because setType reset to 1
      nbColsSub=mat.nbColsSub;  //redo because setType reset to 1

      if (Trace::traceMemory)
        {
          thePrintStream << "LargeMatrix::copy_constructor allocates a new large matrix : " << &values_ << ", "
                         << mat.values_.size() << " non zeros coefficients " << dimValues();
          if (mat.storage_p != 0) thePrintStream << ", storage " << mat.storage_p->name();
          thePrintStream << eol << std::flush;
        }

      values_ = mat.values_;
      storage_p = mat.storage_p;
      if (storage_p == 0) return *this;
      storage_p->objectPlus();
      factorization_=mat.factorization_;
      #ifdef XLIFEPP_WITH_UMFPACK
      umfLU_p=0;
      if (mat.umfLU_p != 0)
       {
        warning("free_warning","LargeMatrix copy, does not copy the umfpackLU object, redo factorisation if required");
        factorization_=_noFactorization;
       }
      #endif
    }
  return *this;
}

template <typename T>
template <typename K>
LargeMatrix<T>::LargeMatrix(const LargeMatrix<K>& mat, bool storageCopy)
  : valueType_(mat.valueType_), strucType(mat.strucType), nbRows(mat.nbRows), nbCols(mat.nbCols), sym(mat.sym),
    nbRowsSub(mat.nbRowsSub), nbColsSub(mat.nbColsSub), name("shallow copy of" + mat.name)
{
  setType(T());
  nbRowsSub=mat.nbRowsSub;  //redo because setType reset to 1
  nbColsSub=mat.nbColsSub;  //redo because setType reset to 1

  if (Trace::traceMemory)
    {
      thePrintStream << "LargeMatrix::copy_constructor allocates a new large matrix : " << &values_ << ", "
                     << mat.values().size() << " non zeros coefficients " << dimValues();
      if (mat.storagep() != 0) thePrintStream << ", storage " << mat.storagep()->name();
      thePrintStream << eol << std::flush;
    }

  values_.resize(mat.values().size());
  typename std::vector<T>::iterator it = values_.begin();
  typename std::vector<K>::const_iterator itm = mat.values().begin();
  for (; it != values_.end(); it++, itm++) *it = *itm;
  storage_p = mat.storagep();
  if (storage_p == 0) return;
  if (storageCopy) storage_p=mat.storagep()->clone();               //storage copy
  storage_p->objectPlus();
  factorization_=mat.factorization_;
  #ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
  if (mat.umfLU_p != 0)
    {
        warning("free_warning","LargeMatrix copy, does not copy the umfpackLU object, redo factorisation if required");
        factorization_=_noFactorization;
    }
  #endif
}

template <typename T>
template <typename K>
LargeMatrix<T>& LargeMatrix<T>::operator=(const LargeMatrix<K>& mat)
{
  if (this != &mat)
    {
      clearall();
      valueType_=mat.valueType_;
      strucType=mat.strucType;
      nbRows=mat.nbRows;
      nbCols=mat.nbCols;
      sym=mat.sym;
      nbRowsSub=mat.nbRowsSub;
      nbColsSub=mat.nbColsSub;
      name="shallow copy of" + mat.name;
      setType(T());
      nbRowsSub=mat.nbRowsSub;  //redo because setType reset to 1
      nbColsSub=mat.nbColsSub;   //redo because setType reset to 1

      if (Trace::traceMemory)
        {
          thePrintStream << "LargeMatrix::copy_constructor allocates a new large matrix : " << &values_ << ", "
                         << mat.values().size() << " non zeros coefficients " << dimValues();
          if (mat.storagep() != 0) thePrintStream << ", storage " << mat.storagep()->name();
          thePrintStream << eol << std::flush;
        }

      values_.resize(mat.values().size());
      typename std::vector<T>::iterator it = values_.begin();
      typename std::vector<K>::const_iterator itm = mat.values().begin();
      for (; it != values_.end(); it++, itm++) *it = *itm;
      storage_p = mat.storagep();
      if (storage_p == 0) return *this;
      storage_p->objectPlus();
      factorization_=mat.factorization_;
      #ifdef XLIFEPP_WITH_UMFPACK
      umfLU_p=0;
      if (mat.umfLU_p != 0)
      {
        warning("free_warning","LargeMatrix copy, does not copy the umfpackLU object, redo factorisation if required");
        factorization_=_noFactorization;
      }
      #endif
    }
  return *this;
}

//constructor from a given storage and a constant value, sy tells if the matrix has symmetry
template <typename T>
LargeMatrix<T>::LargeMatrix(MatrixStorage* ms, const T& val, SymType sy)
  : valueType_(_none), strucType(_scalar), nbRows(0), nbCols(0), sym(sy), nbRowsSub(1), nbColsSub(1)
{
  setType(val);
  init(ms, val, sy);
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  factorization_=_noFactorization;
}

template <typename T>
LargeMatrix<T>::LargeMatrix(MatrixStorage* ms, SymType sy)
  : valueType_(_none), strucType(_scalar), nbRows(0), nbCols(0), sym(sy), nbRowsSub(1), nbColsSub(1)
{
  T val = T(0);
  setType(val);
  init(ms, val, sy);
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  factorization_=_noFactorization;
}

//use this constructor for matrix of matrix, dims are the dimensions of matrix values
template <typename T>
LargeMatrix<T>::LargeMatrix(MatrixStorage* ms, dimPair dims, SymType sy)
  : valueType_(_none), strucType(_scalar), nbRows(0), nbCols(0), sym(sy), nbRowsSub(1), nbColsSub(1)
{
  T val = T();
  resize(val, dims.first, dims.second);
  setType(val);
  init(ms, val, sy);
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  factorization_=_noFactorization;

}

//construct particular matrix such as constant diagonal matrix
template <typename T>
LargeMatrix<T>::LargeMatrix(SpecialMatrix sp, StorageType st, AccessType at, number_t nbr, number_t nbc, const T v)
{
  sym=_noSymmetry;
  factorization_=_noFactorization;
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  if(sp==_idMatrix)
    {
      if(at==_sym) sym=_symmetric;
      number_t m=std::min(nbr,nbc);
      std::vector<std::vector<number_t> > indices(m);  // vector of column indices for each row
      std::vector<std::vector<number_t> >::iterator iti=indices.begin();
      for(number_t k=1; k<=m; k++, iti++) *iti = std::vector<number_t>(1,k);
      MatrixStorage* diagst = createMatrixStorage(st, at, nbr, nbc, indices,"");
      setType(v);
      init(diagst,v,sym);
    }
  else
    {
      where("LargeMatrix<T>::LargeMatrix(SpecialMatrix, StorageType, AccessType, Number, Number, T)");
      error("special_matrix_unexpected", words("matrix", _idMatrix), words("matrix", sp));
    }
}

// destructor : the storage pointer is deleted if no other object shares it
template <typename T>
LargeMatrix<T>::~LargeMatrix()
{
  clearall();
}

// clear values and storage pointer if not shared
template <typename T>
void LargeMatrix<T>::clearall()
{
#ifdef XLIFEPP_WITH_UMFPACK
  if(umfLU_p!=0) delete umfLU_p;
#endif
  clear();
  if(storage_p != 0)
    {
      if(storage_p->numberOfObjects() > 0) storage_p->objectMinus();
      if(storage_p->numberOfObjects() == 0) delete storage_p;
      storage_p =0;
    }
}

//functions to resize matrix or vector values (perhaps to be defined elsewhere)
template <typename T>
void resize(T& v, dimen_t m, dimen_t n)
{
  return;
}
template <typename T>
void resize(std::vector<T>& v, dimen_t m, dimen_t n)
{
  v.resize(m);
  return;
}
template <typename T>
void resize(Matrix<T>& v, dimen_t m, dimen_t n)
{
  v.resize(m * n);
  return;
}

//initialize type of matrix element
template <typename T>
void LargeMatrix<T>::setType(const T& v)
{
  structPair sp = Value::typeOf(v);
  valueType_ = sp.first;
  strucType = sp.second;
  dimPair ds = dimsOf(v);
  nbRowsSub = ds.first;
  nbColsSub = ds.second;
  if(nbRowsSub>1 || nbColsSub>1) strucType=_matrix;
}

// initialize LargeMatrix from a given storage and a constant value, sy tells if the matrix has symmetry
template <typename T>
void  LargeMatrix<T>::init(MatrixStorage* ms, const T& val, SymType sy)
{
  storage_p = ms;
  if(ms == 0) return;       // void storage, allocation delayed

  nbRows = storage_p->nbOfRows();
  nbCols = storage_p->nbOfColumns();
  number_t s = storage_p->size() + 1;
  if(storage_p->accessType() == _sym && sy == _noSymmetry) s+= storage_p->upperPartSize();   // symmetric storage for non symmetric matrix
  if(Trace::traceMemory)
    {
      thePrintStream<<"LargeMatrix::init allocates a new large matrix : "<<&values_<<", "
                    <<s<<" non zeros coefficients "<<dimValues();
      if(storage_p!=0) thePrintStream<<", storage "<<storage_p->name();
      thePrintStream<<eol<<std::flush;
    }
  values_.resize(s, val);
  values_[0] = 0 * val;       // set first to zero (matrix begins at adress 1!)
  storage_p->objectPlus();
}

// constructor with matrix dimensions, storage type and a constant value
// storage pointer is allocated but for other storages than dense storage is not set
template <typename T>
LargeMatrix<T>::LargeMatrix(number_t nr, number_t nc, StorageType sto, AccessType at, const T& val)
  : valueType_(_none), strucType(_scalar), nbRows(nr), nbCols(nc), sym(_noSymmetry)
{
  factorization_=_noFactorization;
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  setType(val);
  allocateStorage(sto, at, val);
  if(at == _sym && sym == _noSymmetry) //force matrix to be symmetric
    {
      warning("largematrix_forcesymmetry");
      sym = _symmetric;
    }
}

// construct "symmetric" matrix from dimensions, storage type, symmetry type and a constant value
// acces type is _sym in that case
template <typename T>
LargeMatrix<T>::LargeMatrix(number_t nr, number_t nc, StorageType sto, SymType sy, const T& val)
  : valueType_(_none), strucType(_scalar), nbRows(nr), nbCols(nc), sym(sy)
{
  factorization_=_noFactorization;
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  setType(val);
  allocateStorage(sto, _sym, val);
  // for skewsymetric matrix, force diagonal to be zero
  if(sy == _skewSymmetric)
    {
      T zero = 0 * val;
      typename std::vector<T>::iterator itb = values_.begin() + 1, it;
      for(it = itb; it < itb + nbRows; it++)
        {
          *it = zero;
        }
    }
  // for skewadjoint matrix, force diagonal to have zero real part
  if(sy == _skewAdjoint)
    {
      T img = (val - conj(val)) / 2.;
      typename std::vector<T>::iterator itb = values_.begin() + 1, it;
      for(it = itb; it < itb + nbRows; it++)
        {
          *it = img;
        }
    }
}

/* construct matrix from file and storage type; matrix in file is either
   in dense format (fst=_dense) :
    A11 A12  ...                            Re(A11) Im(A11) Re(A12) IM(A12) ....
    A21 A22  ...    or if complex matrix    Re(A21) Im(A21) Re(A22) IM(A22) ....
    ...                                     ...
   or in coordinate sparse storage format (fst=_coo) :
    i j Aij         or if complex values    i j Re(Aij) Im(Aij)
   does not work for matrix of matrices!
*/
template <typename T>
LargeMatrix<T>::LargeMatrix(const string_t& fn, StorageType fst, number_t nr, number_t nc,
                            StorageType sto, AccessType at, number_t nsr, number_t nsc)
  : valueType_(_none), strucType(_scalar), nbRows(nr), nbCols(nc), sym(_noSymmetry), nbRowsSub(nsr), nbColsSub(nsc)
{
  factorization_=_noFactorization;
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  setType(T());
  nbRowsSub = nsr;
  nbColsSub = nsc;
  allocateStorage(sto, at);
  if(nsr > 1 || nsc > 1)
    {
      error("largematrix_noloadmatrixofmatrix");
    }
  loadFromFile(fn, fst);
}

// construct "symmetric" matrix from file and storage type
template <typename T>
LargeMatrix<T>::LargeMatrix(const string_t& fn, StorageType fst, number_t nr, StorageType sto,
                            SymType sy, number_t nsr)
  : valueType_(_none), strucType(_scalar), nbRows(nr), nbCols(nr), sym(sy), nbRowsSub(nsr), nbColsSub(nsr)
{
  factorization_=_noFactorization;
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=0;
#endif
  setType(T());
  nbRowsSub = nsr;
  nbColsSub = nsr;
  if(nsr > 1)
    {
      error("largematrix_noloadmatrixofmatrix");
    }
  allocateStorage(sto, _sym);
  loadFromFile(fn, fst);
}

// allocate storage pointer by creating storage object according to storage type and storage access
// if possible (dense storage), the values_ vector is sized
template <typename T>
void LargeMatrix<T>::allocateStorage(StorageType sto, AccessType at, const T& val)
{
  switch(sto)
    {
    case _dense   :
      switch(at)
        {
        case _dual :
          storage_p = new DualDenseStorage(nbRows, nbCols);
          break;
        case _row  :
          storage_p = new RowDenseStorage(nbRows, nbCols);
          break;
        case _col  :
          storage_p = new ColDenseStorage(nbRows, nbCols);
          break;
        case _sym  :
          storage_p = new SymDenseStorage(nbRows);
          break;
        default    :
          error("storage_bad_access", words("access type",at), words("storage type",_dense));
        }
      if(Trace::traceMemory)
        {
          thePrintStream<<"LargeMatrix::allocateStorage allocates a new large matrix : "<<&values_<<", "
                        <<storage_p->size() + 1<<" non zeros coefficients "<<dimValues();
          if(storage_p!=0) thePrintStream<<", storage "<<storage_p->name();
          thePrintStream<<eol<<std::flush;
        }
      values_.resize(storage_p->size() + 1, val); // allocate values vector
      break;
    case _skyline :
      switch(at)
        {
        case _dual :
          storage_p = new DualSkylineStorage(nbRows, nbCols);
          break;
        case _sym  :
          storage_p = new SymSkylineStorage(nbRows);
          break;
        case _row  :
        case _col  :
        default    :
          error("storage_bad_access", words("access type",at), words("storage type",_skyline));
        }
      break;
    case _cs      :
      switch(at)
        {
        case _dual :
          storage_p = new DualCsStorage(nbRows, nbCols);
          break;
        case _row  :
          storage_p = new RowCsStorage(nbRows, nbCols);
          break;
        case _col  :
          storage_p = new ColCsStorage(nbRows, nbCols);
          break;
        case _sym  :
          storage_p = new SymCsStorage(nbRows);
          break;
        default    :
          error("storage_bad_access", words("access type",at), words("storage type",_cs));
        }
      break;
    default :
      where("LargeMatrix<T>::allocateStorage");
      error("storage_not_handled", words("storage type",sto), words("access type",at));
    }
  storage_p->objectPlus();
}

/*---------------------------------------------------------------------------------------------------
  LargeMatrix<T> various utilities
---------------------------------------------------------------------------------------------------*/



//convert matrix storage to skyline storage
template <typename T>
void LargeMatrix<T>::toSkyline()
{
  if(storage_p == 0) error("matrix_nostorage");
  StorageType st = storage_p->storageType();
  if(st == _skyline) return; //matrix is already stored as skyline

  trace_p->push("LargeMatrix<T>::toSkyline");

  //find if storage is already defined
  AccessType at = _dual;
  if(sym != _noSymmetry) at = _sym;
  MatrixStorage*  skysto = 0;
  // MatrixStorage*  skysto = findMatrixStorage(storage_p->stringId, _skyline, at);
  // if(skysto != 0 && (nbRows!=skysto->nbOfRows() || nbCols!=skysto->nbOfColumns())) skysto=0;   //not the same sizes
  // if(skysto == 0)  //construct newskyline storage  ### deactivate by Eric, because re-using existing storage may be hazardous
  // {
  std::vector<number_t> skyrowpointer(storage_p->skylineRowPointer());
  if(at == _sym) skysto = static_cast< MatrixStorage* >(new SymSkylineStorage(skyrowpointer, storage_p->stringId));
  else
    {
      std::vector<number_t> skycolpointer(storage_p->skylineColPointer());
      skysto = static_cast< MatrixStorage* >(new DualSkylineStorage(skyrowpointer, skycolpointer, storage_p->stringId));
    }
  // }
  // skysto->visual(thePrintStream);

  //update values (Note : test memory in future)
  std::vector<T> oldvalues = values_;
  values_.assign(skysto->size() + 1, 0.*oldvalues[0]);
  storage_p->fillSkylineValues(oldvalues, values_, sym); //does not require new skyline storage !!!
  storage_p->objectMinus();
  //Eric's note : we enforce the destruction of the original storage_p object if it is no longer used by other matrix
  if(storage_p->numberOfObjects()==0) delete storage_p;
  storage_p = skysto;
  storage_p->objectPlus();
  trace_p->pop();
}

// Extract and convert matrix storage to UmfPack storage (similar to CSC storage, differ only in the offset of values)
// Only can this function convert from CS storage to UmfPack
template <typename T>
void LargeMatrix<T>::extract2UmfPack(std::vector<int_t>& colPointer, std::vector<int_t>& rowIndex, std::vector<T>& matA) const
{
  trace_p->push("LargeMatrix<T>::extract2UmfPack()");
  if(storage_p==0) error("matrix_nostorage");
  storage_p->toUmfPack(values_, colPointer, rowIndex, matA, sym);
  trace_p->pop();
}

// Get cs col storage pointer and pointer to first value of matrix
//   only available for LargeMatrix  stored as cs col, move before
//   return false if not a cs col largematrix
template <typename T>
bool LargeMatrix<T>::getCsColUmfPack(std::vector<int_t>& colPointer, std::vector<int_t>& rowIndex, const T*& matA) const
{
  trace_p->push("LargeMatrix<T>::getCsColUmfPack");
  if(storage_p==0) error("matrix_nostorage");
  if(storage_p->storageType()!=_cs || storage_p->accessType()!=_col)
    {
      trace_p->pop();  //use extract2Umfpack
      return false;
    }

  ColCsStorage* csto = reinterpret_cast<ColCsStorage*>(storage_p);  //downcast to ColCsStorage
  //move indices from number_t to int_t
  if(csto->size() > number_t(theIntMax))  //check if the last index goes over the int max
    error("conversion_overflow","Number", "Int", theIntMax);
  colPointer.resize(csto->colPointer().size());
  std::vector<number_t>::iterator itn = csto->colPointer().begin();
  std::vector<int_t>::iterator iti=colPointer.begin();
  for(; itn!=csto->colPointer().end(); ++itn, ++iti) *iti = int_t(*itn);
  rowIndex.resize(csto->rowIndex().size());
  itn=csto->rowIndex().begin();
  iti=rowIndex.begin();
  for(; itn!=csto->rowIndex().end(); ++itn, ++iti) *iti = int_t(*itn);

  //get first address of matrix with 1-shift to be compliant with umfpack which starts at 0
  matA = &(values_.at(1));
  trace_p->pop();
  return true;
}

//convert matrix storage to matrix storage ms, general algorithm to be improved
template <typename T>
void LargeMatrix<T>::toStorage(MatrixStorage* ms)
{
  trace_p->push("LargeMatrix<T>::toStorage");
  if(storage_p == 0) error("matrix_nostorage");
  if(ms == 0) error("matrix_nostorage");
  if(storage_p==ms)
    {
      trace_p->pop();  //same storage, nothing to do
      return;
    }
  if(storage_p->accessType()!=_sym && ms->accessType()==_sym) error("nonsym_to_sym", "conversion");

  std::vector<T> oldvalues = values_;        //copy values
  if(Trace::traceMemory)
    {
      thePrintStream<<"LargeMatrix::toStorage re-allocates a large matrix : "<<&values_<<", "
                    <<ms->size()<<" non zero coefficients "<<dimValues();
      if(storage_p!=0) thePrintStream<<", storage "<<storage_p->name();
      thePrintStream<<eol<<std::flush;
    }
  try
    {
      values_.assign(ms->size() + 1, 0.*oldvalues[0]);      //reallocate values_ vector regarding new storage
    }
  catch(std::exception& e)
    {
      std::cout << e.what() << " in LargeMatrix<T>::toStorage : "<<std::endl;
      std::cout << "try to allocate a vector of size "<< ms->size() + 1;
      std::cout << " requiring "<< ((ms->size() + 1)/(1024*1024))*sizeof(T)<<" Mo";
      abort();
    }

  std::vector<number_t> row(1), cols(nbCols), pos, newpos;
  for(number_t c=0; c<nbCols; c++) cols[c]=c+1;
  for(number_t r=1; r<=nbRows; r++)            //travel matrix rows
    {
      row[0]=r;
      storage_p->positions(row, cols, pos, false, sym);
      ms->positions(row, cols, newpos, false, sym);
      std::vector<number_t>::iterator itp=pos.begin(),itnp=newpos.begin();
      if(sym==_noSymmetry || (storage_p->accessType()!=_sym && ms->accessType()!=_sym))  //matrix has no symmetry property or storages are both non symmetric        {
        {
          for(; itp!=pos.end(); itp++, itnp++)           //travel columns
            if(*itp!=0 && *itnp!=0) values_[*itnp]=oldvalues[*itp];
        }
      else //matrix with symmetry
        {
          //copy lower part
          for(number_t c=1; c<=r; itp++, itnp++, c++)           //travel columns
            if(*itp!=0 && *itnp!=0) values_[*itnp]=oldvalues[*itp];
          //copy upper part
          if(storage_p->accessType()==_sym && ms->accessType()!=_sym)
            switch(sym)
              {
              default :
              case _symmetric :
              {
                for(; itp!=pos.end(); itp++, itnp++)           //travel columns
                  if(*itp!=0 && *itnp!=0) values_[*itnp]=oldvalues[*itp];
                break;
              }
              case _selfAdjoint :
              {
                for(; itp!=pos.end(); itp++, itnp++)           //travel columns
                  if(*itp!=0 && *itnp!=0) values_[*itnp]=conj(oldvalues[*itp]);
                break;
              }
              case _skewSymmetric :
              {
                for(; itp!=pos.end(); itp++, itnp++)           //travel columns
                  if(*itp!=0 && *itnp!=0) values_[*itnp]=-oldvalues[*itp];
                break;
              }
              case _skewAdjoint :
              {
                for(; itp!=pos.end(); itp++, itnp++)           //travel columns
                  if(*itp!=0 && *itnp!=0) values_[*itnp]=-conj(oldvalues[*itp]);
                break;
              }
              }
        }
    }
  storage_p->objectMinus();
  if(storage_p->numberOfObjects()<=0) delete storage_p;
  ms->objectPlus();
  storage_p=ms;
  trace_p->pop();
}

// change coefficients to their conjugates
template <typename T>
LargeMatrix<T>& LargeMatrix<T>::toConj()
{
  if(valueType_==_complex)
    {
      typename std::vector<T>::iterator itv=values_.begin();
      for(; itv!=values_.end(); itv++)  *itv = conj(*itv);
    }
  return *this;
}

// create new scalar matrix entry from non scalar matrix entry
template <typename T>
template <typename K>
LargeMatrix<K>* LargeMatrix<T>::toScalar(K bid)
{
  trace_p->push("LargeMatrix<T>::toScalar");
  LargeMatrix<K>* mat= 0;
  if(strucType == _scalar) error("vector_or_matrix");
  MatrixStorage* nsto = storage_p->toScalar(nbRowsSub,nbColsSub);  //create scalar storage from storage
  mat= new LargeMatrix<K>(nsto, sym);  //create matrix
  toScalarEntries(values_,mat->values(),*nsto);  //copy values
  trace_p->pop();
  return mat;
}

/*! copy matrix values to scalar values, generic algorithm using getrow or getcol
  K is the valueType of both matrices, sto is the matrixstorage of the scalar matrix
*/
template <typename T>
template <typename K>
void LargeMatrix<T>::toScalarEntries(const std::vector<Matrix<K> >& val, std::vector<K>& sval, const MatrixStorage& scalsto)
{
  trace_p->push("LargeMatrix<T>::toScalarEntries");

  AccessType acty=accessType();
  number_t rmax = nbRows, cmax = nbCols, srmax = nbRowsSub*nbRows, scmax=nbColsSub*nbCols;

  if(acty !=_col)  //case _row (all matrix),  _sym and _dual (lower part)
    {
      for(number_t r=0; r<nbRows; r++)
        {
          if(acty !=_row) cmax=r+1;
          std::vector<std::pair<number_t, number_t> > colind=getRow(r+1, 1, cmax);
          std::vector<std::pair<number_t, number_t> >::const_iterator itc;
          for(dimen_t rs=1; rs<=nbRowsSub; rs++)
            {
              number_t rss=r*nbRowsSub+rs;
              if(acty !=_row) scmax=rss;
              std::vector<std::pair<number_t, number_t> > scalcolind = scalsto.getRow(sym,rss,1,scmax);
              std::vector<std::pair<number_t, number_t> >::iterator itsc=scalcolind.begin();
              for(itc=colind.begin(); itc!=colind.end(); itc++)
                for(dimen_t cs=1; cs<=nbColsSub && itsc<scalcolind.end(); cs++, itsc++)
                  sval[itsc->second] = val[itc->second](rs,cs);
            }
        }
    }

  if(acty ==_col || acty ==_dual || (acty==_sym && sym==_noSymmetry))  //case _col (all matrix) and _dual (strict upper part part)
    {
      for(number_t c=0; c<nbCols; c++)
        {
          if(acty != _col) rmax=c+1;
          std::vector<std::pair<number_t, number_t> > rowind=getCol(c+1, 1, rmax);
          std::vector<std::pair<number_t, number_t> >::const_iterator itr;
          for(dimen_t cs=1; cs<=nbColsSub; cs++)
            {
              number_t css=c*nbColsSub+cs;
              if(acty !=_col) srmax=css;
              std::vector<std::pair<number_t, number_t> > scalrowind = scalsto.getCol(sym, css, 1, srmax);
              std::vector<std::pair<number_t, number_t> >::iterator itsr=scalrowind.begin();
              for(itr=rowind.begin(); itr!=rowind.end(); itr++)
                for(dimen_t rs=1; rs<=nbRowsSub && itsr<scalrowind.end(); rs++, itsr++)
                  sval[itsr->second] = val[itr->second](rs,cs);
            }
        }
    }
  trace_p->pop();
}

/*! convert current matrix to unsymmetric representation if it has a symmetric one
  convert only matrix using a symmetric access storage,
    if at = _sym  symmetric storage is not changed  (default behaviour)
    if at = _dual we moved to a dual storage
  values vector is extended to store upper part
*/
template <typename T>
void LargeMatrix<T>::toUnsymmetric(AccessType at)
{
  if(storage_p->accessType() != _sym || sym==_noSymmetry) return;  //no conversion
  if(at!=_sym && at!=_dual)
    error("access_unexpected", words("access type",_sym)+" "+words("or")+" "+words("access type", _dual), words("access type", at));
  if(at == _dual)
    {
      MatrixStorage* ms = storage_p->toDual();    //new dual storage from sym storage
      storage_p->objectMinus();
      if(storage_p->numberOfObjects()<=0) delete storage_p;
      ms->objectPlus();
      storage_p = ms;
    }
  //extend values_ vector size
  number_t ls=storage_p->lowerPartSize()+ storage_p->diagonalSize(), as=ls+storage_p->lowerPartSize();
  if(Trace::traceMemory)
    {
      thePrintStream<<"LargeMatrix::toUnsymmetric allocates a new large matrix : "<<&values_<<", "
                    <<as + 1<<" non zeros coefficients "<<dimValues();
      if(storage_p!=0) thePrintStream<<", storage "<<storage_p->name();
      thePrintStream<<eol<<std::flush;
    }
  values_.resize(as+1);
  //"copy" lower part to upper part
  typename std::vector<T>::iterator itlow=values_.begin()+storage_p->diagonalSize()+1,
                                    itup=values_.begin()+ls+1;
  switch(sym)
    {
    case _symmetric     :
      for(; itup!=values_.end(); itup++, itlow++) *itup = *itlow;
      break;
    case _selfAdjoint   :
      for(; itup!=values_.end(); itup++, itlow++) *itup = conj(*itlow);
      break;
    case _skewSymmetric :
      for(; itup!=values_.end(); itup++, itlow++) *itup = - *itlow;
      break;
    case _skewAdjoint   :
      for(; itup!=values_.end(); itup++, itlow++) *itup = -conj(*itlow);
      break;
    default :
      error("not_handled","LargeMatrix<T>::toUnsymmetric()");
    }
  sym=_noSymmetry;
}

// test if the matrix is in fact a diagonal one (true if all non diagonal coefficients vanish)
template<typename T>
bool LargeMatrix<T>::isDiagonal() const
{
  typename std::vector<T>::const_iterator itv1=values_.begin(), itv2, itv;
  if(storage_p->accessType()==_sym || storage_p->accessType()==_dual)
    {
      for(itv=itv1+storage_p->diagonalSize()+1; itv!=values_.end(); itv++)
        if(xlifepp::norm2(*itv)!=0) return false;
      return true;
    }

  std::vector<number_t> diagpos=storage_p->getDiag(); //get diagonal positions
  std::vector<number_t>::iterator itd=diagpos.begin();
  for(; itd!=diagpos.end(); itd++)
    {
      itv2=values_.begin()+*itd;
      for(itv=itv1+1; itv<itv2; itv++)
        if(xlifepp::norm2(*itv)!=0) return false; //a non diagonal coefficient is not equal to zero
      itv1=itv2;
    }
  itv=itv1+1;
  if(itv != values_.end())
    {
      for(; itv!=values_.end(); itv++)
        if(xlifepp::norm2(*itv)!=0) return false;    //a non diagonal coefficient is not equal to zero
    }
  return true;
}

// test if the matrix is in fact the identity (true if all non diagonal coefficients vanish and diagonal coefficients = 1)
template<typename T>
bool LargeMatrix<T>::isId(const double & tol) const
{
  if(!isDiagonal()) return false;
  std::vector<number_t> diagpos=storage_p->getDiag(); //get diagonal positions
  std::vector<number_t>::iterator it=diagpos.begin();
  for(; it!=diagpos.end(); ++it)
    if(xlifepp::norm2(values_[*it]-1.)>tol) return false;
  return true;
}

//round to 0 all coefficients smaller (in norm) than a tolerance (aszero), storage is not modified
template<typename T>
void LargeMatrix<T>::roundToZero(real_t aszero)
{
  typename std::vector<T>::iterator itv=values_.begin()+1;
  for(; itv!=values_.end(); itv++)
    if(xlifepp::norm2(*itv)<aszero) *itv=T(0);
}

/*! assign to current LargeMatrix a LargeMatrix mat (submatrix) given by its ranks (rows, cols) in current matrix
  a void rows or cols  means that the numbering of mat is the same as current matrix
  in that case, if the storage of mat and current matrix are the same, direct summation is used
  else more expansive method is used : use positions function
  !!! target storage has to be update  and T = K has to be available NO CHECK!!!
*/
template <typename T>
template <typename K>
LargeMatrix<T>& LargeMatrix<T>::assign(const LargeMatrix<K>& mat, const std::vector<number_t>& rows, const std::vector<number_t>& cols)
{
  if(strucType != mat.strucType) error("largematrix_mismatch_structure");
  if(dimValues() != mat.dimValues()) error("largematrix_mismatch_size");
  if(sym != _noSymmetry && mat.sym == _noSymmetry) error("nonsym_to_sym", "assign");
  if(sym != _noSymmetry && mat.sym != _noSymmetry && sym != mat.sym) error("largematrix_diff_sym");

  //simple case : same storage and same symmetry
  if(rows.size() == 0 && cols.size() == 0 && storage_p == mat.storagep())
    {
      typename std::vector<T>::iterator itv = values_.begin();
      typename std::vector<K>::const_iterator itm = mat.values().begin();
      if(sym == mat.sym) // same symmetry property, assign all entries
        {
          for(; itv != values_.end(); itv++, itm++) *itv = *itm;
          return *this;
        }
    }

  //general case
  //   use getrow for lower parts and getcol for upper parts when mat is sym or dual access
  //   use getrow when mat is row access ang getcol when mat is col access
  //   no analysis of storage access of target matrix, in some cases it may be slow (for instance col storage in row storage)
  number_t shr=0, shc=0, nbc=mat.nbCols, nbr=mat.nbRows, m, c, r, p;
  if(rows.size()==1) shr=rows[0];
  if(cols.size()==1) shc=cols[0];
  AccessType msat=mat.accessType();
  std::vector<std::pair<number_t, number_t> >::iterator itr, its, itse;
  std::map<number_t,number_t>::iterator itv;
  const std::vector<K>& valmat=mat.values();

  if(msat==_row || msat==_dual || msat==_sym)  //row in row, lower part of mat
    {
      m = nbc;
      for(number_t k=1; k<=nbr; ++k)
        {
          if(rows.size()>1) r=rows[k-1];
          else r=k+shr;
          if(msat!=_row) m=std::min(k, nbc);   //restrict to lower part
          std::vector<std::pair<number_t, number_t> > coladrs=mat.getRow(k,1,m);    // retry col indices of row k in mat storage
          p=coladrs.size();
          if(p>0)
            {
              std::map<number_t,number_t> colmap;   //reorder column indices : index->adress
              if(cols.size()>1)  //implicit col mapping
                for(itr=coladrs.begin(); itr!=coladrs.end(); ++itr) colmap[cols[itr->first-1]]=itr->second;
              else //explicit col mapping
                for(itr=coladrs.begin(); itr!=coladrs.end(); ++itr) colmap[itr->first+shc]=itr->second;
              number_t minc=colmap.begin()->first, maxc=colmap.rbegin()->first;
              std::vector<std::pair<number_t, number_t> > coladrsG=getRow(r,minc,maxc);
              its=coladrsG.begin();
              itse=coladrsG.end();
              for(itv=colmap.begin(); itv!=colmap.end(); ++itv) //assuming valmat is in the same order as values_
                {
                  c=itv->first;
                  while(its->first != c && its != itse) ++its;
                  if(its != itse && its->first == c) values_[its->second] = valmat[itv->second];
                }
            }
        }
    }
  //if(sym == _noSymmetry && (msat==_col || msat==_dual || msat==_sym))   //col in row, strict upper part of mat
  if(accessType() != _sym && (msat==_col || msat==_dual || msat==_sym))   //col in row, strict upper part of mat
    {
      m = nbr;
      number_t d=2;
      if(msat==_col) d=1;
      for(number_t k=d; k<=nbc; ++k)
        {
          if(cols.size()>1) c=cols[k-1];
          else c=k+shc;
          if(msat!=_col) m=std::min(k-1, nbr);   //restrict to strict upper part
          std::vector<std::pair<number_t, number_t> > rowadrs=mat.getCol(k,1,m);    // retry col indices of row k in mat storage
          p=rowadrs.size();
          if(p>0)
            {
              std::map<number_t,number_t> rowmap;
              if(rows.size()>1)  //implicit row mapping
                for(itr=rowadrs.begin(); itr!=rowadrs.end(); ++itr) rowmap[rows[itr->first-1]]=itr->second;
              else //explicit row mapping
                for(itr=rowadrs.begin(); itr!=rowadrs.end(); ++itr) rowmap[itr->first+shr]=itr->second;
              number_t minr=rowmap.begin()->first, maxr=rowmap.rbegin()->first;
              std::vector<std::pair<number_t, number_t> > rowadrsG= getCol(c,minr,maxr);
              its=rowadrsG.begin();
              itse=rowadrsG.end();
              switch(mat.sym)
                {
                case _selfAdjoint :
                {
                  for(itv=rowmap.begin(); itv!=rowmap.end(); ++itv) //assuming valmat is in the same order as values_
                    {
                      r=itv->first;
                      while(its->first != r && its != itse) ++its;
                      if(its != itse && its->first == r) values_[its->second] = conj(valmat[itv->second]);
                    }
                  break;
                }
                case _skewAdjoint:
                {
                  for(itv=rowmap.begin(); itv!=rowmap.end(); ++itv) //assuming valmat is in the same order as values_
                    {
                      r=itv->first;
                      while(its->first != r && its != itse) ++its;
                      if(its != itse && its->first == r) values_[its->second] = -conj(valmat[itv->second]);
                    }
                  break;
                }
                case _skewSymmetric :
                {
                  for(itv=rowmap.begin(); itv!=rowmap.end(); ++itv) //assuming valmat is in the same order as values_
                    {
                      r=itv->first;
                      while(its->first != r && its != itse) ++its;
                      if(its != itse && its->first == r) values_[its->second] = -valmat[itv->second];
                    }
                  break;
                }
                default : //symmetric or non symmetric
                {
                  for(itv=rowmap.begin(); itv!=rowmap.end(); ++itv) //assuming valmat is in the same order as values_
                    {
                      r=itv->first;
                      while(its->first != r && its != itse) ++its;
                      if(its != itse &&  its->first == r) values_[its->second] = valmat[itv->second];
                    }
                }
                }
            }
        }
    }
  return *this;
}

/*! extract sub matrix of current LargeMatrix and return it as a LargeMatrix
    rowIndex : row indices to be extracted (starts from 1)
    colIndex : col indices to be extracted (starts from 1)
    return a pointer to extracted LargeMatrix in the same storage type of current LargeMatrix
*/
template <typename T>
LargeMatrix<T>* LargeMatrix<T>::extract(const std::vector<number_t>& rowIndex, const std::vector<number_t>& colIndex) const
{
  //create new storage
  MatrixStorage* sto=storage_p->extract(rowIndex, colIndex);
  LargeMatrix<T>* lmsub=new LargeMatrix<T>(sto,sym);
  //fill in lmsub using general non optimal algorithm
  number_t nbr=rowIndex.size(), nbc=colIndex.size();
  Vector<number_t> row = trivialNumbering<number_t>(1, nbr), col=trivialNumbering<number_t>(1, nbc);
  std::vector<number_t> adrsub, adrmat;
  sto->positions(row, col, adrsub);
  storage_p->positions(rowIndex, colIndex, adrmat);
  std::vector<number_t>::iterator its=adrsub.begin(), itm=adrmat.begin();
  for(; its!=adrsub.end(); ++its, ++itm) lmsub->values_[*its] = values_[*itm];
  return lmsub;
}

/*! extract block matrix [r1,r2]x[c1,c2] of current LargeMatrix and return it as a LargeMatrix
    return a pointer to extracted LargeMatrix in the same storage type of current LargeMatrix
    r1<=r2,c1<=c2 indices starting from 1
*/
template <typename T>
LargeMatrix<T>* LargeMatrix<T>::extract(number_t r1, number_t r2, number_t c1, number_t c2) const
{
  if (r1 == 0)
  {
    where("LargeMatrix::extract(...)");
    error("is_null", "r1");
  }
  if (r2 < r1)
  {
    where("LargeMatrix::extract(...)");
    error("is_lesser", r2, r1);
  }
  if (r2 > storage_p->nbOfRows())
  {
    where("LargeMatrix::extract(...)");
    error("is_greater", r2, storage_p->nbOfRows());
  }
  if (c1 == 0)
  {
    where("LargeMatrix::extract(...)");
    error("is_null", "c1");
  }
  if (c2 < c1)
  {
    where("LargeMatrix::extract(...)");
    error("is_lesser", c2, c1);
  }
  if (c2 > storage_p->nbOfColumns())
  {
    where("LargeMatrix::extract(...)");
    error("is_greater", c2, storage_p->nbOfColumns());
  }

  StorageType st=storage_p->storageType();
  AccessType at=storage_p->accessType();
  if(st==_dense && (at==_row || at==_col))  //fast method for row or col dense
    {
      number_t nbr=r2-r1+1, nbc=c2-c1+1;
      std::vector<std::vector<number_t> > indices;//fake indices vector
      string_t id = storage_p->stringId+"_"+tostring(r1)+"_"+tostring(r2)+"_"+tostring(c1)+"_"+tostring(c2);
      MatrixStorage* sto=createMatrixStorage(st,at, nbr, nbc, indices,id);
      LargeMatrix<T>* lmsub(sto,sym);
      //fill in lmsub matrix
      if(at==_row)
        {
          typename std::vector<T>::iterator its=lmsub->values_.begin()+1;
          typename std::vector<T>::const_iterator itm;
          for(number_t r=r1; r<=r2; ++r)
            {
              itm=values_.begin()+(r-1)*nbCols+c1;
              for(number_t c=c1; c<=c2; ++c, ++its, ++itm) *its=*itm;
            }
          return lmsub;
        }
      if(at==_col)
        {
          typename std::vector<T>::iterator its=lmsub->values_.begin()+1;
          typename std::vector<T>::const_iterator itm;
          for(number_t c=c1; c<=c2; ++c)
            {
              itm=values_.begin()+(c-1)*nbRows+r1;
              for(number_t r=r1; r<=r2; ++r, ++its, ++itm) *its=*itm;
            }
          return lmsub;
        }
    }
  // call general method
  std::vector<number_t> rowIndex=trivialNumbering(r1,r2), colIndex=trivialNumbering(c1,c2);
  return extract(rowIndex, colIndex);
}


/*---------------------------------------------------------------------------------------------------
  LargeMatrix<T> addition utilities
---------------------------------------------------------------------------------------------------*/
//add LargeMatrix of same type with same storage
template<typename T>
LargeMatrix<T>& LargeMatrix<T>::operator+=(const LargeMatrix<T>& mat)
{
  //if(storage_p != mat.storage_p) error("largematrix_mismatch_storage");
  if(!sameStorage(*storage_p,*mat.storage_p)) error("largematrix_mismatch_storage");
  typename std::vector<T>::iterator it = values_.begin() + 1;
  typename std::vector<T>::const_iterator itm = mat.values_.begin() + 1;
  if(storage_p->accessType() != _sym || sym == mat.sym) //add values
    {
      for(; it != values_.end(); it++, itm++) *it += *itm;
      return *this;
    }
  //case of symmetric storage but the symmetry properties are different, the result has no symmetry
  typename std::vector<T>::iterator itu, ite;
  if(sym != _noSymmetry) //values_ vector has to be extended
    {
      number_t s = values_.size() + storage_p->lowerPartSize();
      if(Trace::traceMemory)
        {
          thePrintStream<<"LargeMatrix::+= re-allocates a large matrix : "<<&values_<<", "
                        <<s<<" non zeros coefficients "<<dimValues();
          if(storage_p!=0) thePrintStream<<", storage "<<storage_p->name();
          thePrintStream<<eol<<std::flush;
        }
      values_.resize(s);           //extend values to lowerPartSize
      it = values_.begin() + 1;
      itu = values_.begin() + 1 + storage_p->diagonalSize() + storage_p->lowerPartSize();
      ite = itu;
      switch(sym)  //assign lower part to upper part
        {
        case _selfAdjoint :
          for(; it != ite; it++, itu++) *itu = conj(*it);
          break;
        case _skewSymmetric :
          for(; it != ite; it++, itu++) *itu = - *it;
          break;
        case _skewAdjoint :
          for(; it != ite; it++, itu++) *itu = - conj(*it);
          break;
        case _symmetric :
        default         :
          for(; it != ite; it++, itu++) *itu = *it;
        }
      sym = _noSymmetry;
    }
  it = values_.begin() + 1;
  ite = values_.begin() + 1 + storage_p->diagonalSize() + storage_p->lowerPartSize();
  itm = mat.values_.begin() + 1;
  for(; it != ite; it++, itm++) *it += *itm; //add lower part
  if(mat.sym == _noSymmetry) //add upper part
    {
      for(; it != values_.end(); it++, itm++) *it += *itm;
      return *this;
    }
  itm = mat.values_.begin() + 1;
  switch(mat.sym)  //add upper part using lower part
    {
    case _selfAdjoint :
      for(; it != values_.end(); it++, itm++) *it += conj(*itm);
      break;
    case _skewSymmetric :
      for(; it != values_.end(); it++, itm++) *it += - *itm;
      break;
    case _skewAdjoint :
      for(; it != values_.end(); it++, itm++) *it += -conj(*itm);
      break;
    case _symmetric :
    default         :
      for(; it != values_.end(); it++, itm++) *it += *itm;
    }
  return *this;
}

template<typename T>
LargeMatrix<T>& LargeMatrix<T>::operator-=(const LargeMatrix<T>& mat)
{
  //if(storage_p != mat.storage_p) error("largematrix_mismatch_storage");
  if(!sameStorage(*storage_p,*mat.storage_p)) error("largematrix_mismatch_storage");
  typename std::vector<T>::iterator it = values_.begin() + 1;
  typename std::vector<T>::const_iterator itm = mat.values_.begin() + 1;
  if(storage_p->accessType() != _sym || sym == mat.sym) //add values
    {
      for(; it != values_.end(); it++, itm++) *it -= *itm;
      return *this;
    }
  //case of symmetric storage but the symmetry properties are different, the result has no symmetry
  typename std::vector<T>::iterator itu, ite;
  if(sym != _noSymmetry) //values_ vector has to be extended
    {
      number_t s = values_.size() + storage_p->lowerPartSize();
      if(Trace::traceMemory)
        {
          thePrintStream<<"LargeMatrix::-= re-allocates a large matrix : "<<&values_<<", "
                        <<s<<" non zeros coefficients "<<dimValues();
          if(storage_p!=0) thePrintStream<<", storage "<<storage_p->name();
          thePrintStream<<eol<<std::flush;
        }
      values_.resize(s);          //extend values to lowerPartSize
      it = values_.begin() + 1;
      itu = values_.begin() + 1 + storage_p->diagonalSize() + storage_p->lowerPartSize();
      ite = itu;
      switch(sym)  //assign lower part to upper part
        {
        case _selfAdjoint :
          for(; it != ite; it++, itu++) *itu = conj(*it);
          break;
        case _skewSymmetric :
          for(; it != ite; it++, itu++) *itu = - *it;
          break;
        case _skewAdjoint :
          for(; it != ite; it++, itu++) *itu = - conj(*it);
          break;
        case _symmetric :
        default         :
          for(; it != ite; it++, itu++) *itu = *it;
        }
      sym = _noSymmetry;
    }
  it = values_.begin() + 1;
  ite = values_.begin() + 1 + storage_p->diagonalSize() + storage_p->lowerPartSize();
  itm = mat.values_.begin() + 1;
  for(; it != ite; it++, itm++) *it -= *itm; //add lower part
  if(mat.sym == _noSymmetry) //add upper part
    {
      for(; it != values_.end(); it++, itm++) *it -= *itm;
      return *this;
    }

  itm = mat.values_.begin() + 1;
  switch(mat.sym)  //add upper part using lower part
    {
    case _selfAdjoint :
      for(; it != values_.end(); it++, itm++) *it -= conj(*itm);
      break;
    case _skewSymmetric :
      for(; it != values_.end(); it++, itm++) *it -= - *itm;
      break;
    case _skewAdjoint :
      for(; it != values_.end(); it++, itm++) *it -= -conj(*itm);
      break;
    case _symmetric :
    default         :
      for(; it != values_.end(); it++, itm++) *it -= *itm;
    }
  return *this;
}

//add to current LargeMatrix some values given by iterators itb to ite
template <typename T>
template<typename iterator>
LargeMatrix<T>& LargeMatrix<T>::add(iterator& itb, iterator& ite)
{
  typename std::vector<T>::iterator itv=values_.begin();
  number_t i=0;
  for(iterator it=itb; itb!=ite && i<values_.size(); it++, itv++, i++) *itv+=*it;
}

// add a matrix (mat) at ranks (rows, cols), values of matrix are stored by row
// when the current matrix has a symetry, only lower part of mat is added
template <typename T>
LargeMatrix<T>& LargeMatrix<T>::add(const std::vector<T>& mat, const std::vector<number_t>& rows, const std::vector<number_t>& cols)
{
  std::vector<number_t> pos;
  storage_p->positions(rows, cols, pos, true, sym); //get adresses of submatrix elements
  typename std::vector<T>::const_iterator itm = mat.begin();
  typename std::vector<T>::iterator itv = values_.begin();
  std::vector<number_t>::iterator itp;
  if(sym == _noSymmetry) //add full matrix mat
    for(itp = pos.begin(); itp != pos.end(); itp++, itm++)
      {
        *(itv + *itp) += *itm;
      }
  else                 //add lower part of mat in "symmetric" matrix
    {
      itp = pos.begin();
      for(number_t i = 0; i < rows.size(); i++, itm += cols.size())
        for(number_t j = 0; j < cols.size(); j++)
          if(rows[i] >= cols[j])
            {
              *(itv + *itp) += *(itm + j);
              itp++;
            }
    }
  return *this;
}

// add a constant to submatrix given by its ranks (rows, cols)
template <typename T>
LargeMatrix<T>& LargeMatrix<T>::add(const T& c, const std::vector<number_t>& rows, const std::vector<number_t>& cols)
{
  std::vector<number_t> pos;
  storage_p->positions(rows, cols, pos, true, sym); //get adresses of submatrix elements (only lower part if current matrix has symmetry)
  typename std::vector<T>::const_iterator itm;
  typename std::vector<T>::iterator itv = values_.begin();
  std::vector<number_t>::iterator itp;
  for(itp = pos.begin(); itp != pos.end(); itp++, itm++)
    {
      *(itv + *itp) += c;
    }
  return *this;
}

/*! add to current LargeMatrix a scaled LargeMatrix mat (submatrix) given by its ranks (rows, cols) in current matrix and a scaling factor (scale)
  a void rows or cols  means that the numbering of mat is the same as current matrix
  in that case, if the storage of mat and current matrix are the same, direct summation is used
  else access to position is used (expansive method)
   no storage update !!!
*/
template <typename T>
template <typename K, typename C>
LargeMatrix<T>& LargeMatrix<T>::add(const LargeMatrix<K>& mat, const std::vector<number_t>& rows,
                                    const std::vector<number_t>& cols, C scale)
{
  if(strucType != mat.strucType) error("largematrix_mismatch_structure");
  if(dimValues() != mat.dimValues()) error("largematrix_mismatch_size");
  if(sym != _noSymmetry && mat.sym == _noSymmetry) error("nonsym_to_sym", "addition");
  if(sym != _noSymmetry && mat.sym != _noSymmetry && sym != mat.sym)
    error("largematrix_diff_sym");

  //simple case : same storage and same symmetry
  if(rows.size() == 0 && cols.size() == 0 && storage_p == mat.storagep())
    {
      typename std::vector<T>::iterator itv = values_.begin();
      typename std::vector<K>::const_iterator itm = mat.values().begin();
      if(sym == mat.sym) // same symmetry property, add all entries
        {
          for(; itv != values_.end(); itv++, itm++) *itv += scale** itm;
          return *this;
        }
    }

  //general case, use positions (no optimisation, to be improved)
  std::vector<number_t> num_r = rows;
  std::vector<number_t> num_c = cols;
  std::vector<number_t>::iterator it;
  //if(num_r.size() == 0) //if rows is void, take all the rows
  if(num_r.size() <= 1) //take all the rows and shift by num_r[0] if exists, no shift else
    {
      number_t n = 1, shr=0;
      if(num_r.size()==1) shr=num_r[0];
      //num_r.resize(nbRows);
      num_r.resize(mat.nbRows);
      n+=shr;
      for(it = num_r.begin(); it != num_r.end(); it++, n++) *it = n;
    }
  //if(num_c.size() == 0) //if cols is void, take all the cols
  if(num_c.size() <= 1) //take all the rows and shift by num_c[0] if exists, no shift else
    {
      number_t n = 1, shc=0;
      if(num_c.size()==1) shc=num_c[0];
      //num_c.resize(nbCols);
      num_c.resize(mat.nbCols);
      n+=shc;
      for(it = num_c.begin(); it != num_c.end(); it++, n++) *it = n;
    }

  //loop on submatrix rows
  number_t r = 1;
  for(it = num_r.begin(); it != num_r.end(); it++, r++)
    {
      //get adresses of row r of matrix mat
      std::vector<number_t> mcols(num_c.size()), mrows(1, r), mpos, cpos;
      std::vector<number_t>::iterator itc, itp, itq;
      number_t c = 1;
      for(itc = mcols.begin(); itc != mcols.end(); itc++, c++) *itc = c;
      mat.storagep()->positions(mrows, mcols, mpos, false, mat.sym);     //get adresses of mat coefficients on row n (zero coefficient have a zero adress)

      //translate row/columns indices into numbering of current matrix and get positions into
      itc = mcols.begin();
      itq = mpos.begin();
      std::vector<number_t>::iterator itn = num_c.begin();
      number_t n = 0;
      for(itp = mpos.begin(); itp != mpos.end(); itp++, itn++)  //keep non zero addresses
        {
          if(*itp != 0)
            {
              *itc = *itn;
              *itq = *itp;
              itc++;
              itq++;
              n++;
            }
        }
      mcols.resize(n);
      mpos.resize(n);
      mrows[0] = *it;
      positions(mrows, mcols, cpos, false); //get addresses of row n in current matrix

      // do combination
      itp = cpos.begin(), itq = mpos.begin();
      if(sym == _noSymmetry)
        {
          switch(mat.sym)
            {
            case _skewSymmetric :
              for(itc = mcols.begin(); itc != mcols.end() ; itc++, itp++, itq++)
                if(*it >= *itc) values_[*itp] += scale * mat(*itq); //lower triangular part
                else values_[*itp] -= scale * mat(*itq);            //upper triangular part
              break;
            case _selfAdjoint :
              for(itc = mcols.begin(); itc != mcols.end() ; itc++, itp++, itq++)
                if(*it >= *itc) values_[*itp] += scale * mat(*itq); //lower triangular part
                else values_[*itp] += scale * conj(mat(*itq));       //upper triangular part
              break;
            case _skewAdjoint :
              for(itc = mcols.begin(); itc != mcols.end() ; itc++, itp++, itq++)
                if(*it >= *itc) values_[*itp] += scale * mat(*itq); //lower triangular part
                else values_[*itp] -= scale * conj(mat(*itq));       //upper triangular part
              break;
            default : //no symmetry or symmetric
              for(itc = mcols.begin(); itc != mcols.end(); itc++, itp++, itq++)
                values_[*itp] += scale * mat(*itq);
            }
        }
      else if(mat.sym == sym) //same symmetry property
        {
          for(itc = mcols.begin(); itc != mcols.end(); itc++, itp++, itq++)
            if(*it >= *itc) values_[*itp] += scale * mat(*itq);
        }
      else
        error("largematrix_mismatch_sym");
    }
  return *this;
}

// combine cols of matrix : c2 = c2 + a *c1 (indices start from 1) with update storage option
// use getCols
template <typename T>
void LargeMatrix<T>::addColToCol(number_t c1, number_t c2, complex_t a, bool updateStorage)
{
  ScalarType b = complexToT<ScalarType>(a);
  std::set<number_t> rowsc1= storage_p->getRows(c1), rowsc2=storage_p->getRows(c2);
  rowsc1.insert(rowsc2.begin(),rowsc2.end());
  if(rowsc1.size()==rowsc2.size())  // same set of rows (compatible storage) do combination
    {
      std::vector<std::pair<number_t, number_t> > adrs1=getCol(c1), adrs2=getCol(c2);
      std::vector<std::pair<number_t, number_t> >::iterator it1=adrs1.begin(), it2=adrs2.begin();
      std::map<number_t,number_t> mapadrs2;
      for(; it2!=adrs2.end(); it2++)  mapadrs2[it2->first] = it2->second;
      for(; it1!=adrs1.end(); it1++)  values_[mapadrs2[it1->first]] += b *values_[it1->second];
      return;
    }
  if(!updateStorage)
    {
      where("LargeMatrix<T>::addColToCol(...)");
      error("storage_not_updated");
    }
  //update storage (expansive  operation)
  error("not_yet_implemented",string_t("LargeMatrix<T>::addColToCol(...)")+string_t(" ")+words("with dynamic update storage"));
}

// combine rows of matrix : r2 = r2 + a *r1 (indices start from 1) with update storage option
// use getRows
template <typename T>
void LargeMatrix<T>::addRowToRow(number_t r1, number_t r2, complex_t a, bool updateStorage)
{
  ScalarType b = complexToT<ScalarType>(a);
  std::set<number_t> colsr1= storage_p->getCols(r1), colsr2=storage_p->getCols(r2);
  colsr1.insert(colsr2.begin(),colsr2.end());
  if(colsr1.size()==colsr2.size())  // compatible storage, do combination
    {
      std::vector<std::pair<number_t, number_t> > adrs1=getRow(r1), adrs2=getRow(r2);
      std::vector<std::pair<number_t, number_t> >::iterator it1=adrs1.begin(), it2=adrs2.begin();
      std::map<number_t,number_t> mapadrs2;
      for(; it2!=adrs2.end(); it2++)  mapadrs2[it2->first] = it2->second;
      for(; it1!=adrs1.end(); it1++)  values_[mapadrs2[it1->first]] += b *values_[it1->second];
      return;
    }
  if(!updateStorage)
    {
      where("LargeMatrix<T>::addColToCol(...)");
      error("storage_not_updated");
    }
  //update storage (expansive  operation)
  error("not_yet_implemented",string_t("LargeMatrix<T>::addRowToRow(...)")+string_t(" ")+words("with dynamic update storage"));
}

// norm of a LargeMatrix (Frobenius norm)
template <typename T>
real_t LargeMatrix<T>::norm2() const
{
  real_t n=0.;
  typename std::vector<T>::const_iterator itm=values_.begin()+1;
  for(; itm!=values_.end(); ++itm)
    {
      real_t t = xlifepp::norm2(*itm);
      n+=t*t;
    }
  if(storage_p->accessType()==_sym) //add the strict lower part
    {
      itm=values_.begin()+1+storage_p->diagonalSize();
      for(; itm!=values_.end(); ++itm)
        {
          real_t t = xlifepp::norm2(*itm);
          n+=t*t;
        }
    }
  return std::sqrt(n);
}

template <typename T>
inline real_t norm2(const LargeMatrix<T>& A)
{
  return A.norm2();
}

// norm infinite of a LargeMatrix
template <typename T>
real_t LargeMatrix<T>::norminfty() const
{
  real_t n=0.;
  typename std::vector<T>::const_iterator itm=values_.begin()+1;
  for(; itm!=values_.end(); ++itm)
    {
      real_t t = xlifepp::norminfty(*itm);
      if(t>n) n=t;
    }
  return n;
}

template <typename T>
inline real_t norminfty(const LargeMatrix<T>& A)
{
  return A.norminfty();
}

// partial norm of column c, only components in range [r1,r2]
template <typename T>
real_t LargeMatrix<T>::partialNormOfCol(number_t c, number_t r1, number_t r2) const
{
  real_t n=0.;
  std::vector<std::pair<number_t, number_t> > rowadrs=getCol(c,r1,r2);
  std::vector<std::pair<number_t, number_t> >::iterator itra=rowadrs.begin();
  for(; itra!=rowadrs.end(); itra++) n += xlifepp::norm2(values_[itra->second]);
  return std::sqrt(n);
}

//! get col c (>=1) as vector
template <typename T>
std::vector<T> LargeMatrix<T>::col(number_t c) const
{
  std::vector<T> vcol(nbRows,T());
  std::vector<std::pair<number_t, number_t> > rowadr=getCol(c); // get (row indices, adresses) of col c
  std::vector<std::pair<number_t, number_t> >::iterator ita=rowadr.begin();
  for(; ita!=rowadr.end(); ++ita) vcol[ita->first-1]=values_[ita->second];
  return vcol;
}

//! get row r (>=1) as vector
template <typename T>
std::vector<T> LargeMatrix<T>::row(number_t r) const
{
  std::vector<T> vrow(nbCols,T());
  std::vector<std::pair<number_t, number_t> > coladr=getRow(r); // get (row indices, adresses) of col c
  std::vector<std::pair<number_t, number_t> >::iterator ita=coladr.begin();
  for(; ita!=coladr.end(); ++ita) vrow[ita->first-1]=values_[ita->second];
  return vrow;
}

// delete rows r1,..., r2 (may be expansive for some storage)
template <typename T>
void LargeMatrix<T>::deleteRows(number_t r1, number_t r2)
{
  if(storage_p!=0)
    {
      storage_p->deleteRows(r1, r2, values_,sym);
      nbRows=storage_p->nbOfRows();
    }
}
// delete cols c1,..., c2 (may be expansive for some storage)
template <typename T>
void LargeMatrix<T>::deleteCols(number_t c1, number_t c2)
{
  if(storage_p!=0)
    {
      storage_p->deleteCols(c1, c2, values_,sym);
      nbCols=storage_p->nbOfColumns();
    }
}

// set to 0 cols/rows cr1,..., cr2 (index start from 1), no col/row deletion
// generic algorithm use getCol/getRow of storage,
// may be slightly improved by using storage child versions, not yet done
template <typename T>
void LargeMatrix<T>:: setColToZero(number_t c1, number_t c2)
{
  if(c1==0)
    {
      c1=1;
      c2=nbCols;
    }
  if(c2==0)
    {
      c2=nbCols;
    }
  if(c1>c2) return;
  std::vector<std::pair<number_t, number_t> > adrs;
  std::vector<std::pair<number_t, number_t> >::iterator ita;
  for(number_t c=c1; c<=c2; c++)
    {
      adrs = getCol(c, 1, 0);  //get adresses of all non zero of column c
      if(sym ==_noSymmetry)
        for(ita=adrs.begin(); ita!=adrs.end(); ita++) values_[ita->second]=T(0);
      else //matrix with symmetry property
        for(ita=adrs.begin(); ita!=adrs.end(); ita++)
          if(ita->first >= c) values_[ita->second]=T(0);
    }
}

template <typename T>
void LargeMatrix<T>:: setRowToZero(number_t r1, number_t r2)
{
  if(r1==0)
    {
      r1=1;
      r2=nbRows;
    }
  if(r2==0)
    {
      r2=nbRows;
    }
  if(r1>r2) return;
  std::vector<std::pair<number_t, number_t> > adrs;
  std::vector<std::pair<number_t, number_t> >::iterator ita;
  for(number_t r=r1; r<=r2; r++)
    {
      adrs = getRow(r, 1, 0);  //get adresses of all non zero of row c
      if(sym ==_noSymmetry)
        for(ita=adrs.begin(); ita!=adrs.end(); ita++) values_[ita->second]=T(0);
      else //matrix with symmetry property
        for(ita=adrs.begin(); ita!=adrs.end(); ita++)
          if(ita->first <= r) values_[ita->second]=T(0);
    }
}

/*! copy in current matrix values of mat located at mat_adrs at position cur_adrs
  it is assumed that mat_adrs and cur_adrs have the same size (address map!)
*/
template <typename T>
void LargeMatrix<T>::copyVal(const LargeMatrix<T>& mat, const std::vector<number_t>& mat_adrs, const std::vector<number_t>& cur_adrs)
{
  if (cur_adrs.size() != mat_adrs.size())
  {
    where("LargeMatrix<T>::copyVal(...)");
    error("bad_size", "address map", mat_adrs.size(), cur_adrs.size());
  }
  std::vector<number_t>::const_iterator itma=mat_adrs.begin(), itca=cur_adrs.begin();
  number_t sma=mat.values_.size()-1, sca=values_.size()-1;
  for (; itma!=mat_adrs.end(); itma++, itca++)
  {
    if (*itma > sma)
    {
      where("LargeMatrix<T>::copyVal(...)");
      error("index_out_of_range", "itma", 0, sma);
    }
    if (*itca >sca)
    {
      where("LargeMatrix<T>::copyVal(...)");
      error("index_out_of_range", "itca", 0, sca);
    }
    values_[*itca]=mat.values_[*itma];
  }
}

/*---------------------------------------------------------------------------------------------------
  LargeMatrix<T> print utilities
---------------------------------------------------------------------------------------------------*/
template <typename T>
void LargeMatrix<T>::print(std::ostream& os) const
{
  if(theVerboseLevel == 0)
    {
      return;
    }
  number_t nnz = nbNonZero();
  string_t vt = words("value", valueType_) + " " + words("structure", strucType);
  string_t sto = "undefined storage";
  if(storage_p != 0)
    {
      sto = words("access type", storage_p->accessType()) + " " + words("storage type", storage_p->storageType());
    }
  os << message("largematrix_header", words("symmetry", sym), vt, nbRows, nbCols, sto, nnz)<<eol;
  if(theVerboseLevel > 0 && storage_p != 0)
    {
      storage_p->printEntries(os, values_, theVerboseLevel, sym);
      number_t n=rowPermutation_.size();
      if(n>0)
        {
          os<<"row permutation : [";
          if(n<=2*theVerboseLevel)
            {
              for(number_t i=0; i<n; i++) os<<rowPermutation_[i]<<" ";
            }
          else
            {
              for(number_t i=0; i<theVerboseLevel; i++) os<<rowPermutation_[i]<<" ";
              os<<"... ";
              for(number_t i=n-theVerboseLevel; i<n; i++) os<<rowPermutation_[i]<<" ";
            }
          os<<"]"<<eol;
        }
      n=colPermutation_.size();
      if(n>0)
        {
          os<<"column permutation : [";
          if(n<=2*theVerboseLevel)
            {
              for(number_t i=0; i<n; i++) os<<colPermutation_[i]<<" ";
            }
          else
            {
              for(number_t i=0; i<theVerboseLevel; i++) os<<colPermutation_[i]<<" ";
              os<<"... ";
              for(number_t i=n-theVerboseLevel; i<n; i++) os<<colPermutation_[i]<<" ";
            }
          os<<"]"<<eol;
        }
    }
}

//! output stream operator
template <typename T>
std::ostream& operator<<(std::ostream& os, const LargeMatrix<T>& mat)
{
  mat.print(os);
  return os;
}

//visualize storage on output stream (like spy on matlab)
template <typename T>
void LargeMatrix<T>::viewStorage(std::ostream& os) const
{
  storage_p->visual(os);
}

/*---------------------------------------------------------------------------------------------------
  save to file utilities
---------------------------------------------------------------------------------------------------*/
/* save matrix to a file along different format
   only _dense and _coo storage type are allowed
   the file format being consistent with Matlab load functions
   if encodeName is true, the name of file includes some additionnal matrix informations
    newname = name(m_n_SorageType_ValueType_StructType_p_q).ext
   (m,n) size of matrix, (p,q) size of submatrices if StrucType=_matrix
*/

template <typename T>
void LargeMatrix<T>::saveToFile(const string_t& fn, StorageType st, bool encodeName) const
{
  string_t fen = fn;
  if(encodeName)
    {
      fen = encodeFileName(fn, st);
    }
  std::ofstream of(fen.c_str());
  if(of.fail())
    {
      error("file_failopen", "LargeMatrix<T>::saveToFile", fen);
    }
  of.precision(fullPrec);
  switch(st)
    {
    case _dense :
      storage_p->printDenseMatrix(of, values_, sym);
      break;
    case _coo   :
      storage_p->printCooMatrix(of, values_, sym);
      break;
    default :
      error("largematrix_nosavematrix", words("storage type", st));
    }
  of.close();
}

template <typename T>
void LargeMatrix<T>::saveToFile(const char* fn, StorageType st, bool encodeName) const
{
  string_t fen = fn;
  saveToFile(fen,st,encodeName);
}

//encode file name with additionnal matrix informations
//  n=name.ext -> name(m_n_SorageType_ValueType_StructType_p_q).ext
// (m,n) size of matrix, (p,q) size of submatrices if StrucType=_matrix
template <typename T>
string_t  LargeMatrix<T>::encodeFileName(const string_t& fn, StorageType st) const
{
  number_t pos = fn.find_last_of('.');
  string_t fr, ext = "", stn = "dense", vt = "real";
  if(valueType_ == _complex)
    {
      vt = "complex";
    }
  if(st == _coo)
    {
      stn = "coo";
    }
  if(pos != string_t::npos)
    {
      fr = fn.substr(0, pos);
      ext = fn.substr(pos);
    }
  else
    {
      fr = fn;
    }
  fr += "(" + tostring(nbRows) + "_" + tostring(nbCols) + "_" + stn + "_" + vt;
  if(strucType == _scalar)
    {
      fr += "_scalar)" + ext;
    }
  else
    {
      fr += "_scalar_" + tostring(nbRowsSub) + "_" + tostring(nbColsSub) + ")" + ext;
    }
  return fr;
}

/*---------------------------------------------------------------------------------------------------
 load matrix from a file either
   in dense format (s=_dense)     :   A11 A12 ...  A1n       re(A11) im(A11)  re(A21) im(A21) ...
                                      A21 A22 ...  A2n       re(A21) im(A21)  re(A22) im(A22) ...
                    ....                   ...
   in coordinate format (st=_coo) :   i j Aij (one per row) or if complex values : i j re(Aij) im(Aij)
 only space separator
 storage type, acces type, value type, struct type and symmetry of the matrix have to be already set
 calls the virtual storage_p->loadFromFileDense or storage_p->loadFromFileCoo functions
---------------------------------------------------------------------------------------------------*/
template <typename T>
void LargeMatrix<T>::loadFromFile(const string_t& fn, StorageType stf)
{
  trace_p->push("LargeMatrix::loadFromFile");
  std::ifstream ifs(fn.c_str());
  if(ifs.fail()) error("file_failopen", "LargeMatrix::loadFromFile", fn);

  //check compatibility by analysing first line
  string_t line;
  std::getline(ifs, line);
  line = trim(line);
  std::stringstream ss(line);
  real_t r;
  number_t nbc = 0;
  while(!ss.eof())
    {
      ss >> r;
      nbc++;
    }
  number_t nbcw = 0, nbcwr = 0;
  if(stf == _dense)
    {
      nbcw = nbCols * nbColsSub;
      nbcwr = 2 * nbcw;
    }
  else
    {
      nbcw = 3;
      nbcwr = 4;
    }
  if((valueType_ == _real && nbc != nbcw) || (valueType_ == _complex && nbc != nbcw && nbc != nbcwr))
    {
      ss.str(string_t());
      ss.clear();
      ss << "LargeMatrix::loadFromFile(" << fn << ", " << words("storage type", stf) << ")";
      error("bad_dim", ss.str(), nbcw, nbc);
    }
  bool rAsC = (valueType_ == _complex && nbc == nbcw);
  //call child functions
  ifs.clear();
  ifs.seekg(0, std::ios::beg);  //go to the beginning of file
  storage_p->clear();
  switch(stf)
    {
    case _dense:
      storage_p->loadFromFileDense(ifs, values_, sym, rAsC);
      break;
    case _coo:
      storage_p->loadFromFileCoo(ifs, values_, sym, rAsC);
      break;
    default:
      where("LargeMatrix::loadFromFile");
      error("storage_not_handled", words("storage type", stf), words("access type", _noAccess));
    }
  ifs.close();
  trace_p->pop();
}

/*======================================================================================================
  special matrix product used by svd compression methods
  LargeMatrix * MatCol or MatCol * LargeMatrix or LargeMatrix * MatRow or MatRow * LargeMatrix
  =====================================================================================================*/
/*! do the product R = A * M where
       A is any largeMatrix
       pM pointer to the first element of the dense col matrix
       pR pointer to the first element of the dense col result matrix
       p the number of cols of M
       assuming the number of rows of M is equal to the number of cols of A
*/
template<typename T>
void LargeMatrix<T>::multMatrixCol(T* pM, T* pR, number_t p) const
{
  number_t m=numberOfRows(), n=numberOfCols();
  T* pRk=pR, *pMk=pM;
  typename std::vector<T>::const_iterator itAb=values_.begin()+1, itA;
  if(storageType()==_dense) //fast algorithm
    {
      if(accessType()==_row)
        {
          for(number_t j=0; j<p; ++j)
            for(number_t i=0; i<m; ++i, ++pRk)
              {
                T rij=T();
                itA=itAb+i*n;
                pMk=pM+j*n;
                for(number_t k=0; k<n; ++k, ++itA, ++pMk) rij+=*itA**pMk;
                *pRk=rij;
              }
          return;
        }
      if(accessType()==_col)
        {
          for(number_t j=0; j<p; ++j)
            for(number_t i=0; i<m; ++i, ++pRk)
              {
                T rij=T();
                itA=itAb+i;
                pMk=pM+j*n;
                for(number_t k=0; k<n; ++k, itA+=m, ++pMk) rij+=*itA**pMk;
                *pRk=rij;
              }
          return;
        }
    }
  //general case : use matrix vector product
  for(number_t k=0; k<p; ++k, pRk+=n, pMk+=n)
    multMatrixVector(*this, pMk, pRk);
}

/*! do the product R = M * A where
       A is any largeMatrix
       pM pointer to the first element of the dense col matrix
       pR pointer to the first element of the dense col result matrix
       p the number of rows of M
       assuming the number of rows of M is equal to the number of cols of A
*/
template<typename T>
void LargeMatrix<T>::multLeftMatrixCol(T* pM, T* pR, number_t p) const
{
  number_t m=numberOfRows(), n=numberOfCols();
  T* pRk=pR, *pMk=pM;
  typename std::vector<T>::const_iterator itAb=values_.begin()+1, itA;
  if(storageType()==_dense) //fast algorithm
    {
      if(accessType()==_row)
        {
          for(number_t j=0; j<n; ++j)
            for(number_t i=0; i<p; ++i, ++pRk)
              {
                T rij=T();
                itA=itAb+j;
                pMk=pM+i;
                for(number_t k=0; k<m; ++k, itA+=n, pMk+=p) rij+=*pMk ** itA;
                *pRk=rij;
              }
          return;
        }
      if(accessType()==_col)
        {
          for(number_t j=0; j<n; ++j)
            for(number_t i=0; i<p; ++i, ++pRk)
              {
                T rij=T();
                itA=itAb+j*m;
                pMk=pM+i;
                for(number_t k=0; k<m; ++k, itA++, pMk+=p) rij+=*pMk ** itA;
                *pRk=rij;
              }
          return;
        }
    }
  //general case : use matrix vector product
  std::vector<T> x(n), y(m);
  typename std::vector<T>::iterator itx;
  for(number_t i=0; i<p; ++i)
    {
      pMk=pM+i;
      for(itx=x.begin(); itx!=x.end(); ++itx, pMk+=p) *itx = *pMk;
      multVectorMatrix(*this, x, y);
      pRk=pR+i;
      for(itx=y.begin(); itx!=y.end(); ++itx, pRk+=p) *pRk = *itx;
    }
}


/*! do the product R = A * M where
       A is any largeMatrix
       pM pointer to the first element of the dense row matrix
       pR pointer to the first element of the dense row result matrix
       p the number of cols of M
       assuming the number of rows of M is equal to the number of cols of A
*/
template<typename T>
void LargeMatrix<T>::multMatrixRow(T* pM, T* pR, number_t p) const
{
  number_t m=numberOfRows(), n=numberOfCols();
  T* pRk=pR, *pMk=pM;
  typename std::vector<T>::const_iterator itAb=values_.begin()+1, itA;
  if(storageType()==_dense) //fast algorithm
    {
      if(accessType()==_row)
        {
          for(number_t i=0; i<m; ++i)
            for(number_t j=0; j<p; ++j, ++pRk)
              {
                T rij=T();
                itA=itAb+i*n;
                pMk=pM+j;
                for(number_t k=0; k<n; ++k, ++itA, pMk+=p) rij+=*itA**pMk;
                *pRk=rij;
              }
          return;
        }
      if(accessType()==_col)
        {
          for(number_t i=0; i<m; ++i)
            for(number_t j=0; j<p; ++j, ++pRk)
              {
                T rij=T();
                itA=itAb+i;
                pMk=pM+j;
                for(number_t k=0; k<n; ++k, itA+=m, pMk+=p) rij+=*itA**pMk;
                *pRk=rij;
              }
          return;
        }
    }
  //general case : use matrix vector product
  std::vector<T> x(n), y(m);
  typename std::vector<T>::iterator itx;
  for(number_t j=0; j<p; ++j)
    {
      pMk=pM+j;
      for(itx=x.begin(); itx!=x.end(); ++itx, pMk+=p) *itx = *pMk;
      multMatrixVector(*this, x, y);
      pRk=pR+j;
      for(itx=y.begin(); itx!=y.end(); ++itx, pRk+=p) *pRk = *itx;
    }
}

/*! do the product R = M * A where
       A is any largeMatrix
       pM pointer to the first element of the dense row matrix
       pR pointer to the first element of the dense row result matrix
       p the number of rows of M
       assuming the number of rows of M is equal to the number of cols of A
*/
template<typename T>
void LargeMatrix<T>::multLeftMatrixRow(T* pM, T* pR, number_t p) const
{
  number_t m=numberOfRows(), n=numberOfCols();
  T* pRk=pR, *pMk=pM;
  typename std::vector<T>::const_iterator itAb=values_.begin()+1, itA;
  if(storageType()==_dense) //fast algorithm
    {
      if(accessType()==_row)
        {
          for(number_t i=0; i<p; ++i)
            for(number_t j=0; j<n; ++j, ++pRk)
              {
                T rij=T();
                itA=itAb+j;
                pMk=pM+i*m;
                for(number_t k=0; k<m; ++k, itA+=n, pMk++) rij+=*pMk**itA;
                *pRk=rij;
              }
          return;
        }
      if(accessType()==_col)
        {
          for(number_t i=0; i<p; ++i)
            for(number_t j=0; j<n; ++j, ++pRk)
              {
                T rij=T();
                itA=itAb+j*m;
                pMk=pM+i*m;
                for(number_t k=0; k<m; ++k, itA++, pMk++) rij+=*pMk**itA;
                *pRk=rij;
              }
          return;
        }
    }
  //general case : use matrix vector product
  std::vector<T> x(m), y(n);
  typename std::vector<T>::iterator itx;
  for(number_t i=0; i<p; ++i)
    {
      pMk=pM+i*m;
      for(itx=x.begin(); itx!=x.end(); ++itx, pMk++) *itx = *pMk;
      multVectorMatrix(*this, x, y);
      pRk=pR+i*n;
      for(itx=y.begin(); itx!=y.end(); ++itx, pRk++) *pRk = *itx;
    }
}

/*---------------------------------------------------------------------------------------------------
  extern operations for LargeMatrix
  for each functionnality, consistant template version
  and consistant specializations (real_t/complex_t Scalar/Vector/Matrix) are proposed
---------------------------------------------------------------------------------------------------*/
/*!product with factorized matrix
      no permutation  : L*U*X or L*D*Lt*X or (L*D*L*)*X
      row permutation  (PA=LU) :  inv(P)*L*U*X or inv(P)*L*D*Lt*X or inv(P)*(L*D*L*)*X
      col permutation  (AQ=LU) :  L*U*inv(Q)*X or L*D*Lt*inv(Q)*X or (L*D*L*)*inv(Q)*X
      row and col permutation  (PAQ=LU) :  inv(P)*L*U*inv(Q)*X or inv(P)*L*D*Lt*inv(Q)*X or inv(P)*(L*D*L*)*inv(Q)*X
   called by multMatrixVector functions
*/
template <typename T, typename V, typename R>
void multFactMatrixVector(const LargeMatrix<T>& mat, const std::vector<V>& vec, std::vector<R>& res)
{
  number_t n=vec.size();
  res.resize(n);
  typename std::vector<V>::const_iterator itv=vec.begin();
  typename std::vector<R>::iterator itr;
  std::vector<R> r1(n), r2(n);
  for(itr=r1.begin(); itr!=r1.end(); ++itr, ++itv) *itr=*itv;
  if(mat.colPermutation_.size()>0) permuteInv(r1,r1,mat.colPermutation_);
  switch(mat.factorization_)
    {
    case _ilu:
    case _lu:
      mat.storage_p->upperMatrixVector(mat.values_, r1, r2, mat.sym);
      mat.storage_p->lowerD1MatrixVector(mat.values_, r2, res, mat.sym);
      break;
    case _ldlt:
      mat.storage_p->upperD1MatrixVector(mat.values_, r1, r2, mat.sym);
      mat.storage_p->diagonalMatrixVector(mat.values_, r2, r1, mat.sym);
      mat.storage_p->lowerD1MatrixVector(mat.values_, r1, res, mat.sym);
      break;
    case _ldlstar:
      mat.storage_p->upperD1MatrixVector(mat.values_, r1, r2, mat.sym);
      mat.storage_p->diagonalMatrixVector(mat.values_, r2, r1, mat.sym);
      for(itr=r1.begin(); itr!=r1.end(); ++itr) *itr=conj(*itr); //conjugate res conj(A)*X= conj(A*conj(X))
      mat.storage_p->lowerD1MatrixVector(mat.values_, r1, res, mat.sym);
      for(itr=res.begin(); itr!=res.end(); ++itr) *itr=conj(*itr); //conjugate res
      break;
    default :
      where("multFactMatrixVector(LargeMatrix, vector, vector)");
      error("wrong_factorization_type", words("factorization type", mat.factorization_));
    }
  if(mat.rowPermutation_.size()>0) permuteInv(res,res,mat.rowPermutation_);
}

/*!product with factorized matrix A=LU, or A=LDLt or A =LDL*
      no permutation  : X*L*U = (Ut*Lt*Xt)t
      row permutation  (PA=LU) : X*inv(P)*L*U = Ut*Lt*inv(P)t*Xt =  Ut*Lt*P*Xt
      col permutation  (AQ=LU) : X*L*U*inv(Q) = inv(Q)t*Ut*Lt*Xt =  Q*Ut*Lt*Xt
      row and col permutation  (PAQ=LU) :  X*inv(P)*L*U*inv(Q) = inv(Q)t*Ut*Lt*inv(P)t*Xt =  Q*Ut*Lt*P*Xt
   called by multMatrixVector
*/
template <typename T, typename V, typename R>
void multVectorFactMatrix(const LargeMatrix<T>& mat, const std::vector<V>& vec, std::vector<R>& res)
{

  error("free_error","product Vector * Factorized Matrix is not yet available");

//    ### DOES NOT WORK, lowerVectorMatrix, upperVectorMatrix, ... not yet available n storage class###
//    number_t n=vec.size();
//    res.resize(n);
//    typename std::vector<V>::const_iterator itv=vec.begin();
//    typename std::vector<R>::iterator itr;
//    std::vector<R> r1(n), r2(n);
//    for(itr=r1.begin();itr!=r1.end();++itr, ++itv) *itr=*itv;
//    if(mat.rowPermutation_.size()>0) permute(r1,r1,mat.rowPermutation_);
//    switch(mat.factorization_)
//    {
//    case _lu: mat.storage_p->vectorLowerD1Matrix(mat.values_, r1, r2, mat.sym);
//              mat.storage_p->vectorUpperMatrix(mat.values_, r2, res, mat.sym);
//
//        break;
//    case _ldlt: mat.storage_p->vectorLowerD1Matrix(mat.values_, r1, res, mat.sym);
//                mat.storage_p->diagonalMatrixVector(mat.values_, r2, r1, mat.sym);
//                mat.storage_p->lowerD1MatrixVector(mat.values_, r1, r2, mat.sym);
//
//        break;
//    case _ldlstar: mat.storage_p->vectorLowerD1Matrix(mat.values_, r1, res, mat.sym);
//                   mat.storage_p->diagonalMatrixVector(mat.values_, r2, r1, mat.sym);
//                   for(itr=r1.begin();itr!=r1.end();++itr) *itr=conj(*itr);  //conjugate res conj(A)*X= conj(A*conj(X))
//                   mat.storage_p->lowerD1MatrixVector(mat.values_, r1, r2, mat.sym);
//                   for(itr=res.begin();itr!=res.end();++itr) *itr=conj(*itr);  //conjugate res
//        break;
//    default :  error("free_error","factorization type not handled in multFactMatrixVector");
//    }
//    if(mat.colPermutation_.size()>0) permute(res,res,mat.colPermutation_);
}

//! templated  mat<T> * vec<T>
template <typename T>
void multMatrixVector(const LargeMatrix<T>& mat, const std::vector<T>& vec, std::vector<T>& res)
{
  trace_p->push("multMatrixVector");
  if (vec.size() != mat.nbCols) error("largematrix_mismatch_dim");
  if (res.size() < mat.nbRows) res.resize(mat.nbRows, 0);
  if (mat.factorization_==_noFactorization) mat.storage_p->multMatrixVector(mat.values_, vec, res, mat.sym);
  else multFactMatrixVector(mat, vec, res);
  trace_p->pop();
}

//! templated  mat<T> * vec<V> (pointer form)
template <typename T, typename V, typename R>
void multMatrixVector(const LargeMatrix<T>& mat, V* vp, R* rp)
{
  mat.storage_p->multMatrixVector(mat.values_, vp, rp, mat.sym);
}

//! templated  vec<V> * mat<T>  (pointer form)
template <typename T, typename V, typename R>
void multVectorMatrix(const LargeMatrix<T>& mat, V* vp, R* rp)
{
  mat.storage_p->multVectorMatrix(mat.values_, vp, rp, mat.sym);
}

//! templated  vec<V> * mat<T>  (pointer form)
template <typename T, typename V, typename R>
void multVectorMatrix(V* vp, const LargeMatrix<T>& mat, R* rp)
{
  mat.storage_p->multVectorMatrix(mat.values_, vp, rp, mat.sym);
}

template <typename T>
std::vector<T> operator*(const LargeMatrix<T>& mat, const std::vector<T>& vec)
{
  trace_p->push("LargeMatrix * vector");
  std::vector<T> res(mat.nbRows, T());
  multMatrixVector(mat, vec, res);
  trace_p->pop();
  return res;
}

//! templated  mat<Matrix<T> > * vec<Vector<T> >
template <typename T>
void multMatrixVector(const LargeMatrix<Matrix<T> >& mat, const std::vector<Vector<T> >& vec, std::vector<Vector<T> >& res)
{
  trace_p->push("multMatrixVector(LargeMatrix<Matrix>, vector<Vector>, vector<Vector>)");
  if(vec.size() != mat.nbCols) error("largematrix_mismatch_dim");
  if(res.size() < mat.nbRows)  res.resize(mat.nbRows, Vector<T>(mat.nbRowsSub, T()));
  if(mat.factorization_==_noFactorization) mat.storage_p->multMatrixVector(mat.values_, vec, res, mat.sym);
  else error("largematrix_factorized", mat.name);
  trace_p->pop();
}

template <typename T>
std::vector<Vector<T> > operator*(const LargeMatrix<Matrix<T> >& mat, const std::vector<Vector<T> >& vec)
{
  std::vector<Vector<T> > res(mat.nbRows, Vector<T>(mat.nbRowsSub, T()));
  multMatrixVector(mat, vec, res);
  return res;
}

// specializations mat * vec
std::vector<complex_t> operator*(const LargeMatrix<real_t>& mat, const std::vector<complex_t>& vec);
std::vector<complex_t> operator*(const LargeMatrix<complex_t>& mat, const std::vector<real_t>& vec);
std::vector<Vector<complex_t> > operator*(const LargeMatrix<Matrix<real_t> >& mat, const std::vector<Vector<complex_t> >& vec);
std::vector<Vector<complex_t> > operator*(const LargeMatrix<Matrix<complex_t> >& mat, const std::vector<Vector<real_t> >& vec);

//templated  vec<T> * mat<T>
template <typename T>
void multVectorMatrix(const LargeMatrix<T>& mat, const std::vector<T>& vec, std::vector<T>& res)
{
  trace_p->push("multVectorMatrix");
  if(vec.size() != mat.nbRows) error("largematrix_mismatch_dim");
  if(res.size() < mat.nbCols)  res.resize(mat.nbCols, 0);
  if(mat.factorization_==_noFactorization) mat.storage_p->multVectorMatrix(mat.values_, vec, res, mat.sym);
  else multVectorFactMatrix(mat, vec, res);
  trace_p->pop();
}

//templated  vec<T> * mat<T>
template <typename T>
void multVectorMatrix(const std::vector<T>& vec, const LargeMatrix<T>& mat, std::vector<T>& res)
{
  trace_p->push("multVectorMatrix");
  if(vec.size() != mat.nbRows) error("largematrix_mismatch_dim");
  if(res.size() < mat.nbCols)  res.resize(mat.nbCols, 0);
  if(mat.factorization_==_noFactorization) mat.storage_p->multVectorMatrix(mat.values_, vec, res, mat.sym);
  else multVectorFactMatrix(mat, vec, res);
  trace_p->pop();
}

template <typename T> template <typename K>
LargeMatrix<T>& LargeMatrix<T>::operator*=(const K& v)
{
  typename std::vector<T>::iterator it = values_.begin();
  for(; it != values_.end(); it++)  *it *= v ;
  return *this;
}

template<typename S1, typename S2>
void multInverMatrixVector(const LargeMatrix<S1>& mat, std::vector<S2>& vec,
                           std::vector<typename Conditional<NumTraits<S1>::IsComplex, S1, S2>::type>& res,
                           FactorizationType fac)
{
  switch(fac)
    {
    case _lu:
      mat.luSolve(vec, res);
      break;
    case _ldlt:
      mat.ldltSolve(vec, res);
      break;
    case _ldlstar:
      mat.ldlstarSolve(vec, res);
      break;
    case _umfpack:
      mat.umfluSolve(vec, res);
      break;
    default:
      error("largematrix_not_factorized", mat.name);
    }
}

template <typename T> template <typename K>
LargeMatrix<T>& LargeMatrix<T>::operator/=(const K& v)
{
  typename std::vector<T>::iterator it = values_.begin();
  for(; it != values_.end(); it++)
    {
      *it /= v ;
    }
  return (*this);
}

template <typename T>
std::vector<T> operator*(const std::vector<T>& vec, const LargeMatrix<T>& mat)
{
  std::vector<T> res(mat.nbCols, T());
  multVectorMatrix(mat, vec, res);
  return (res);
}

//templated  vec<Vector<T> > * mat<Matrix<T> >
template <typename T>
void multVectorMatrix(const LargeMatrix<Matrix<T> >& mat, const std::vector<Vector<T> >& vec, std::vector<Vector<T> >& res)
{
  trace_p->push("multVectorMatrix(LargeMatrix<Matrix>, vector<Vector>, vector<Vector>)");
  if(vec.size() != mat.nbRows) error("largematrix_mismatch_dim");
  if(res.size() < mat.nbCols)  res.resize(mat.nbCols, Vector<T>(mat.nbColsSub, T()));
  if(mat.factorization_==_noFactorization) mat.storage_p->multVectorMatrix(mat.values_, vec, res, mat.sym);
  else  error("largematrix_factorized", mat.name);
  trace_p->pop();
}

//templated  vec<Vector<T> > * mat<Matrix<T> >
template <typename T>
void multVectorMatrix(const std::vector<Vector<T> >& vec, const LargeMatrix<Matrix<T> >& mat, std::vector<Vector<T> >& res)
{
  trace_p->push("multVectorMatrix(vector<Vector>, LargeMatrix<Matrix, vector<Vector>)");
  if(vec.size() != mat.nbRows) error("largematrix_mismatch_dim");
  if(res.size() < mat.nbCols)  res.resize(mat.nbCols, Vector<T>(mat.nbColsSub, T()));
  if(mat.factorization_==_noFactorization) mat.storage_p->multVectorMatrix(mat.values_, vec, res, mat.sym);
  else  error("largematrix_factorized", mat.name);
  trace_p->pop();
}

template <typename T>
std::vector<Vector<T> > operator*(const std::vector<Vector<T> >& vec, const LargeMatrix<Matrix<T> >& mat)
{
  std::vector<Vector<T> > res(mat.nbCols, Vector<T>(vec.begin()->size(), T()));
  multVectorMatrix(mat, vec, res);
  return res;
}

// specialize version vec * mat
std::vector<complex_t> operator*(const std::vector<complex_t>& vec, const LargeMatrix<real_t>& mat);
std::vector<complex_t> operator*(const std::vector<real_t>& vec, const LargeMatrix<complex_t>& mat);
std::vector<Vector<complex_t> > operator*(const std::vector<Vector<complex_t> >& vec, const LargeMatrix<Matrix<real_t> >& mat);
std::vector<Vector<complex_t> > operator*(const std::vector<Vector<real_t> >& vec, const LargeMatrix<Matrix<complex_t> >& mat);

// ---------------------------------------------------------------------
//  Factorisation LDLt, LDL*, LU of LargeMatrix
// ---------------------------------------------------------------------
// LDLt Factorization of a symmetric matrix
template <typename T>
void LargeMatrix<T>::ldltFactorize()
{
  trace_p->push("LargeMatrix::ldlt");
  if(_symmetric != sym) storage_p->noFactorization("L.D.Lt");
  StorageType st = storage_p->storageType();
  if(st!=_skyline && st!=_dense) storage_p->noFactorization("L.D.Lt");
  storage_p->ldlt(values_, values_);
  factorization_=_ldlt;
  trace_p->pop();
}

template<typename S>
void ldltFactorize(LargeMatrix<S>& mat)
{
  mat.ldltFactorize();
}

// LDL* Factorization of a Hermitian positive definite matrix
template <typename T>
void LargeMatrix<T>::ldlstarFactorize()
{
  trace_p->push("LargeMatrix::ldlstar");
  if(_selfAdjoint != sym) storage_p->noFactorization("L.D.LStar");
  if(_skyline != (storage_p->storageType())) storage_p->noFactorization("L.D.LStar");
  storage_p->ldlstar(values_, values_);
  factorization_=_ldlstar;
  trace_p->pop();
}

template<typename S>
void ldlstarFactorize(LargeMatrix<S>& mat)
{
  mat.ldlstarFactorize();
}

/*! do LU factorization of skyline or dense matrix
    LU with row permutation is available only for dense matrices
    if matrix is stored as symmetric it is extended to store upper part (toUnsymmetric function)
*/
template<typename T>
void LargeMatrix<T>::luFactorize(bool withPermutation)
{
  trace_p->push("luFactorization");
  if(_skyline == (storage_p->storageType()))
    {
      toUnsymmetric(accessType());
      storage_p->lu(values_, values_,sym);
      factorization_=_lu;
      trace_p->pop();
      return;
    }
  else if(_dense == storage_p->storageType())
    {
      if(_sym != storage_p->accessType()) toUnsymmetric(_dual);
      if(withPermutation) storage_p->lu(values_, values_, rowPermutation_, sym);
      else storage_p->lu(values_, values_, sym);
      factorization_=_lu;
      trace_p->pop();
      return;
    }
  storage_p->noFactorization("L.U");
  trace_p->pop();
}

template<typename S>
void luFactorize(LargeMatrix<S>& mat, bool withPermutation=true)
{
  mat.luFactorize(withPermutation);
}

/*  Factorization incomplete LU without permutation */

template<typename T>
void LargeMatrix<T>::iluFactorize()
{
  trace_p->push("iluFactorization");
  if( _cs != (storage_p->storageType()) ) storage_p->noFactorization("IL.U");
  toUnsymmetric(accessType());
  storage_p->ilu(values_, values_,sym);
  factorization_=_ilu;
  trace_p->pop();
}

template<typename S>
void iluFactorize(LargeMatrix<S>& mat)
{
  mat.iluFactorize();
}


/*  Factorization incomplete LDLt without permutation */
// ---------------------------------------------------------------------
// iLDLt Factorization of a symmetric matrix
template <typename T>
void LargeMatrix<T>::ildltFactorize()
{
  trace_p->push("LargeMatrix::ildlt");
  if(_symmetric != sym) storage_p->noFactorization("iL.D.Lt");
  StorageType st = storage_p->storageType();
  if(st!=_skyline && st!=_dense && st!=_cs) storage_p->noFactorization("iL.D.Lt");
  storage_p->ildlt(values_, values_);
  factorization_=_ildlt;
  trace_p->pop();
}

template<typename S>
void ildltFactorize(LargeMatrix<S>& mat)
{
  mat.ildltFactorize();
}

/*  Factorization incomplete LLt without permutation */
// ---------------------------------------------------------------------
// iLLt Factorization of a symmetric matrix positive definate
template <typename T>
void LargeMatrix<T>::illtFactorize()
{
  trace_p->push("LargeMatrix::illt");
  if(_selfAdjoint != sym && _symmetric != sym) storage_p->noFactorization("iL.Lt");
  StorageType st = storage_p->storageType();
  if(st!=_skyline && st!=_dense && st!=_cs) storage_p->noFactorization("iL.Lt");
  storage_p->illt(values_, values_);
  factorization_=_illt;
  trace_p->pop();
}

template<typename S>
void illtFactorize(LargeMatrix<S>& mat)
{
  mat.illtFactorize();
}

// iLDL* Factorization of a hermitian matrix
template <typename T>
void LargeMatrix<T>::ildlstarFactorize()
{
  trace_p->push("LargeMatrix::ildlstar");
  if( _selfAdjoint != sym && _symmetric != sym ) storage_p->noFactorization("iL.D.LstarR");
  StorageType st = storage_p->storageType();
  if(st!=_skyline && st!=_dense && st!=_cs) storage_p->noFactorization("iL.D.Lstar");
  storage_p->ildlstar(values_, values_);
  factorization_=_ildlstar;
  trace_p->pop();
}

template<typename S>
void ildlstarFactorize(LargeMatrix<S>& mat)
{
  mat.ildlstarFactorize();
}

/*  Factorization incomplete LLstar without permutation */
// ---------------------------------------------------------------------
// iLLt Factorization of a symmetric matrix positive definate
template <typename T>
void LargeMatrix<T>::illstarFactorize()
{
  trace_p->push("LargeMatrix::illstar");
  if( _selfAdjoint != sym && _symmetric != sym) storage_p->noFactorization("iL.Lstarr");
  StorageType st = storage_p->storageType();
  if(st!=_skyline && st!=_dense && st!=_cs) storage_p->noFactorization("iL.Lstar");
  storage_p->illstar(values_, values_);
  factorization_=_illstar;
  trace_p->pop();
}

template<typename S>
void illstarFactorize(LargeMatrix<S>& mat)
{
  mat.illstarFactorize();
}

template<typename T>
void LargeMatrix<T>::umfpackFactorize()
{
  trace_p->push("umfpackFactorization");
#ifdef XLIFEPP_WITH_UMFPACK
  umfLU_p=new UmfPackLU<LargeMatrix<T> >(*this);
  factorization_=_umfpack;
#else //no umfpack
  error("xlifepp_without_umfpack");
#endif // XLIFEPP_WITH_UMFPACK
  trace_p->pop();
}

template<typename S>
void umfpackFactorize(LargeMatrix<S>& mat)
{
  mat.umfpackFactorize();
}

// ---------------------------------------------------------------------
//  Solver for LDLt, LDL*, LU system, matrix has to be factorized before
// ---------------------------------------------------------------------
// Template version of solver with matrix factorization LLT
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::lltSolve(std::vector<S1>& vec1, std::vector<S2>& vec2) const
{
  trace_p->push("LargeMatrix::lltSolve");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::lltSolve", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->sorLowerSolver(values_, vec1, vec2,1.);
  storage_p->sorUpperSolver(values_, vec2, vec2, 1.,sym);
  trace_p->pop();
}
// Template version of solver with matrix factorization LDLT
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::ldltSolve(std::vector<S1>& vec1, std::vector<S2>& vec2) const
{
  trace_p->push("LargeMatrix::ldltSolve");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::ldltSolve", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->lowerD1Solver(values_, vec1, vec2);
  storage_p->diagonalSolver(values_, vec2, vec2);
  storage_p->upperD1Solver(values_, vec2, vec2, sym);
  trace_p->pop();
}

// Template version of solver with matrix factorization LDL*
template <typename T> template<typename S>
void LargeMatrix<T>::ldlstarSolve(std::vector<S>& vec1, std::vector<T>& vec2) const
{
  trace_p->push("LargeMatrix::ldlstarSolve");
  if(sym != _selfAdjoint) storage_p->noSolver("L.D.L*");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::ldlstarSolve", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->lowerD1Solver(values_, vec1, vec2);
  storage_p->diagonalSolver(values_, vec2, vec2);
  storage_p->upperD1Solver(values_, vec2, vec2, sym);
  trace_p->pop();
}

// Template version of solver of LU-factorized matrix
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::luSolve(std::vector<S1>& vec1, std::vector<S2>& vec2) const
{
  trace_p->push("LargeMatrix::luSolve");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::luSolve", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  if(rowPermutation_.size()>0) permute(vec1,vec1,rowPermutation_);
  storage_p->lowerD1Solver(values_, vec1, vec2);
  storage_p->upperSolver(values_, vec2, vec2, sym);
  if(colPermutation_.size()>0) permute(vec2,vec2,colPermutation_);
  trace_p->pop();
}

// Template version of solver of umfpack LU-factorized matrix
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::umfluSolve(std::vector<S1>& vec1, std::vector<S2>& vec2) const
{
#ifdef XLIFEPP_WITH_UMFPACK
  trace_p->push("LargeMatrix::umfluSolver");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::umfluSolve", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  if(umfLU_p==0) error("largematrix_not_factorized",name);
  umfLU_p->solve(vec1,vec2);
  trace_p->pop();
#else
  error("xlifepp_without_umfpack");
#endif
}

// Template version of sor Matrix Vector
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::sorLowerMatrixVector(const std::vector<S1>& vec1, std::vector<S2>& vec2, const real_t w) const
{
  trace_p->push("LargeMatrix::sorLowerMatrixVector");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::sorLowerMatrixVector", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->sorLowerMatrixVector(values_, vec1, vec2, w);
  trace_p->pop();
}

// Template version of sor Matrix Vector
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::sorDiagonalMatrixVector(const std::vector<S1>& vec1, std::vector<S2>& vec2, const real_t w) const
{
  trace_p->push("LargeMatrix::sorDiagonalMatrixVector");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::sorDiagonalMatrixVector", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->sorDiagonalMatrixVector(values_, vec1, vec2, w);
  trace_p->pop();
}

// Template version of sor Matrix Vector
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::sorUpperMatrixVector(const std::vector<S1>& vec1, std::vector<S2>& vec2, const real_t w) const
{
  trace_p->push("LargeMatrix::sorUpperMatrixVector");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::sorUpperMatrixVector", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->sorUpperMatrixVector(values_, vec1, vec2, w, sym);
  trace_p->pop();
}

// Template version of sor solver
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::sorLowerSolve(const std::vector<S1>& vec1, std::vector<S2>& vec2, const real_t w) const
{
  trace_p->push("LargeMatrix::sorLowerSolver");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::sorLowerSolver", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->sorLowerSolver(values_, vec1, vec2, w);
  trace_p->pop();
}

// Template version of sor solver
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::sorDiagonalSolve(const std::vector<S1>& vec1, std::vector<S2>& vec2, const real_t w) const
{
  trace_p->push("LargeMatrix::sorDiagonalSolver");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::sorDiagonalSolver", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->sorDiagonalSolver(values_, vec1, vec2, w);
  trace_p->pop();
}

// Template version of sor solver
template <typename T> template<typename S1, typename S2>
void LargeMatrix<T>::sorUpperSolve(const std::vector<S1>& vec1, std::vector<S2>& vec2, const real_t w) const
{
  trace_p->push("LargeMatrix::sorUpperSolver");
  number_t n=vec1.size();
  if(n != numberOfRows())  error("bad_dim", "LargeMatrix::sorUpperSolver", vec1.size(), vec2.size());
  if(n != vec2.size()) vec2.resize(n);
  storage_p->sorUpperSolver(values_, vec1, vec2, w, sym);
  trace_p->pop();
}

// ---------------------------------------------------------------------
//  interface to UmfPack to solve AX=B
//      current LargeMatrix : A
//      vec : right hand side B
//      res : solution X
//      soa : boolean flag to stop on abnormal situation in UmfPack
//      return the reciproqual condition number rcond = 1/||A||*||inv(A)|| in norm 1
//  Note : these functions call implicitely LU factorisation, see luFactorize for LU factorisation with umfpack
// ---------------------------------------------------------------------
template <typename T> template<typename S>
real_t LargeMatrix<T>::umfpackSolve(const std::vector<S>& vec, std::vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType,
                                    S>::type>& res, bool soa)
{
#ifdef XLIFEPP_WITH_UMFPACK
  UmfPackSolver umPSolver;  //umfPack solver object
  return umPSolver(*this, vec, res, soa);
#else
  error("xlifepp_without_umfpack");
  return real_t(0.); // dummy return
#endif
}

//multiple right hand side, vec has to be of the same type
template <typename T> template<typename S>
real_t LargeMatrix<T>::umfpackSolve(const std::vector<std::vector<S>*>& bs, std::vector<std::vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType,
                                    S>::type>* >& xs, bool soa)
{
#ifdef XLIFEPP_WITH_UMFPACK
  UmfPackSolver umPSolver;  //umfPack solver object
  return umPSolver(*this, bs, xs, soa);
#else
  error("xlifepp_without_umfpack");
  return real_t(0.); // dummy return
#endif
}

// other interface to UmfPack to solve AX=B
//      colPointer, rowIndex : sparse strage information
//      values : values of matrix
//      vec : right hand side B
//      res : solution X
//      soa : boolean flag to stop on abnormal situation in UmfPack
//      return the reciproqual condition number rcond = 1/||A||*||inv(A)|| in norm 1
template <typename T> template<typename S>
real_t LargeMatrix<T>::umfpackSolve(const std::vector<int_t>& colPointer, const std::vector<int_t>& rowIndex, const std::vector<T>& values, const std::vector<S>& vec,
                                    std::vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, S>::type>& res, bool soa)
{
#ifdef XLIFEPP_WITH_UMFPACK
  UmfPackSolver umPSolver; //umfPack solver object
  return umPSolver(nbRows, colPointer, rowIndex, &values[0], vec, res, _A, soa);
#else
  error("xlifepp_without_umfpack");
  return real_t(0.); // dummy return
#endif
}

/* ---------------------------------------------------------------------
   solve triangular systems
       LD1*x = b   (lowerD1Solve)
       U*b = b     (upperSolve)
       x*LD1 = b   (lowerD1LeftSolve)
       x*U = b     (upperLeftSolve)
 ---------------------------------------------------------------------*/
// solve triangular system LD1*x = b
template <typename T>
std::vector<T>& LargeMatrix<T>::lowerD1Solve(std::vector<T>& b, std::vector<T>& x) const
{
  if(nbRows!=b.size())
    {
      where("lowerD1Solve(vector, vector)");
      error("largematrix_mismatch_dim");
    }
  if(b.size() != nbRows) error("bad_dim", "LargeMatrix::lowerD1Solve(vector, vector)", b.size(), nbRows);
  if(nbCols != x.size()) x.resize(nbCols);
  if(rowPermutation_.size()>0) permute(b, b, rowPermutation_);
  storage_p->lowerD1Solver(values_, b, x);
  if(colPermutation_.size()>0) permute(x, x, colPermutation_);
  return x;
}

// solve triangular system U*x = b
template <typename T>
std::vector<T>& LargeMatrix<T>::upperSolve(std::vector<T>& b, std::vector<T>& x) const
{
  if(b.size() != nbRows)  error("bad_dim", "LargeMatrix::upperSolve(vector, vector)", b.size(), nbRows);
  if(nbCols != x.size()) x.resize(nbCols);
  if(rowPermutation_.size()>0) permute(b, b, rowPermutation_);
  storage_p->upperSolver(values_, b, x, sym);
  if(colPermutation_.size()>0) permute(x, x, colPermutation_);
  return x;
}

// solve triangular system x*LD1 = b
template <typename T>
std::vector<T>& LargeMatrix<T>::lowerD1LeftSolve(std::vector<T>& b, std::vector<T>& x) const
{
  if(b.size() != nbCols)  error("bad_dim", "LargeMatrix::lowerD1LeftSolve(vector, vector)", b.size(), nbRows);
  if(nbRows != x.size()) x.resize(nbRows);
  if(colPermutation_.size()>0) permute(b, b, colPermutation_);
  storage_p->lowerD1LeftSolver(values_, b, x);
  if(rowPermutation_.size()>0) permute(x, x, rowPermutation_);
  return x;
}

// solve triangular system x*U = b
template <typename T>
std::vector<T>& LargeMatrix<T>::upperLeftSolve(std::vector<T>& b, std::vector<T>& x) const
{
  if(b.size() != nbCols)  error("bad_dim", "LargeMatrix::upperLeftSolve(vector, vector)", b.size(), nbRows);
  if(nbRows != x.size()) x.resize(nbRows);
  if(colPermutation_.size()>0) permute(b, b, colPermutation_);
  storage_p->upperLeftSolver(values_, b, x, sym);
  if(rowPermutation_.size()>0) permute(x, x, rowPermutation_);
  return x;
}

/* ---------------------------------------------------------------------
   solve triangular systems
       LD1*X = M   (lowerD1Solve)
       U*X = M     (upperSolve)
       X*LD1 = M   (lowerD1LeftSolve)
       X*U = M     (upperLeftSolve)
   the matrix result X is stored as a col dense LargeMatrix
 ---------------------------------------------------------------------*/
// solve triangular system LD1*X = M
template <typename T>
LargeMatrix<T>& LargeMatrix<T>::lowerD1Solve(const LargeMatrix<T>& M, const LargeMatrix<T>& X) const
{
  if(nbRows!=M.nbRows)
    {
      where("lowerD1Solve(LargeMatrix,LargeMatrix)");
      error("largematrix_mismatch_dim");
    }
  X=LargeMatrix(nbCols, M.nbCols, _dense, _col);
  typename std::vector<T>::iterator itX = X.values_.begin()+1, itx;
  std::vector<T> x(nbCols), b(M.nbRows);
  for(number_t k=1; k<=M.nbCols; ++k)
    {
      b=M.col(k);
      lowerD1Solve(b, x,1.);
      itx=x.begin();
      for(number_t i=0; i<nbCols; i++,++itX, ++itx) *itX=*itx;
    }
  return X;
}

template <typename T>
LargeMatrix<T>& LargeMatrix<T>::upperSolve(const LargeMatrix<T>& M, const LargeMatrix<T>& X) const
{
  if(nbRows!=M.nbRows)
    {
      where("lowerD1Solve(LargeMatrix,LargeMatrix)");
      error("largematrix_mismatch_dim");
    }
  X=LargeMatrix(nbCols, M.nbCols, _dense, _col);
  typename std::vector<T>::iterator itX = X.values_.begin()+1, itx;
  std::vector<T> x(nbCols), b(M.nbRows);
  for(number_t k=1; k<=M.nbCols; ++k)
    {
      b=M.col(k);
      upperSolve(b, x,1.);
      itx=x.begin();
      for(number_t i=0; i<nbCols; i++,++itX, ++itx) *itX=*itx;
    }
  return X;
}

template <typename T>
LargeMatrix<T>& LargeMatrix<T>::lowerD1LeftSolve(const LargeMatrix<T>& M, const LargeMatrix<T>& X) const
{
  if(nbCols!=M.nbCols)
    {
      where("lowerD1Solve(LargeMatrix,LargeMatrix)");
      error("largematrix_mismatch_dim");
    }
  X=LargeMatrix(M.nbRows, nbRows, _dense, _col);
  typename std::vector<T>::iterator itX = X.values_.begin()+1, itx;
  std::vector<T> x(M.nbRows), b(M.nbRows);
  for(number_t k=1; k<=M.nbCols; ++k)
    {
      b=M.col(k);
      lowerD1LeftSolve(b, x,1.);
      itx=x.begin();
      for(number_t i=0; i<M.nbRows; i++,++itX, ++itx) *itX=*itx;
    }
  return X;

}

template <typename T>
LargeMatrix<T>& LargeMatrix<T>::upperLeftSolve(const LargeMatrix<T>& M, const LargeMatrix<T>& X) const
{
  if(nbCols!=M.nbCols)
    {
      where("lowerD1Solve(LargeMatrix,LargeMatrix)");
      error("largematrix_mismatch_dim");
    }
  X=LargeMatrix(M.nbRows, nbRows, _dense, _col);
  typename std::vector<T>::iterator itX = X.values_.begin()+1, itx;
  std::vector<T> x(M.nbRows), b(M.nbRows);
  for(number_t k=1; k<=M.nbCols; ++k)
    {
      b=M.col(k);
      upperLeftSolve(b, x,1.);
      itx=x.begin();
      for(number_t i=0; i<M.nbRows; i++,++itX, ++itx) *itX=*itx;
    }
  return X;
}


// ---------------------------------------------------------------------
//  Construct a diagonal Matrix from a size, storage and a specific value
//  st, at : storage type
//  nbr, nbc : number of rows and columns (may be different)
//
//  return diagonal matrix
// ---------------------------------------------------------------------
template <typename T>
LargeMatrix<T> diagonalMatrix(StorageType st, AccessType at, number_t nbr, number_t nbc, const T v)
{
  return LargeMatrix<T>(_diagonal, st, at, nbr, nbc, v);
}

// ---------------------------------------------------------------------
//  Construct a diagonal Matrix from a Matrix with a specific value
//  \param mat prototype matrix
//  \param v value of diagonal
//  \return diagonal matrix
// ---------------------------------------------------------------------
template <typename T>
LargeMatrix<T> diagonalMatrix(const LargeMatrix<T>& mat, const T v)
{
  trace_p->push("diagonalMatrix");
  if(mat.nbRows != mat.nbCols)  error("mat_nonsquare","diagoalMatrix()", mat.nbRows, mat.nbCols);
  if(0 == mat.nbRows)           error("is_void","mat");
  // Same as copy constructor, there is an issue with the number of reference to a storage
  // Only when a largeMatrix is created by loading from a file, does this problem happen.
  // The following is only a temporary solution, the root solution ... not this time!!!
  if(0 == mat.storage_p->numberOfObjects())
    {
      mat.storage_p->objectPlus();
    }
  LargeMatrix<T> diagMat(mat.storage_p, T(), mat.sym);
  diagMat.storage_p->setDiagValue(diagMat.values_, v);
  trace_p->pop();
  return diagMat;
}

// ---------------------------------------------------------------------
//  Construct an identity Matrix from a Matrix
//  \param mat prototype matrix
//  \return identity matrix
// ---------------------------------------------------------------------
template <typename T>
LargeMatrix<T> identityMatrix(const LargeMatrix<T>& mat)
{
  return diagonalMatrix(mat, T(1.0));
}

// ---------------------------------------------------------------------
//  Add two largeMatrix and store result into the third.
//     The two matrices must share the same storage.
//     The result matrix will point to the same storage after the summation
//  \param matA first matrix
//  \param matB second matrix
//  \param matC result matrix which share the same storage.
// ---------------------------------------------------------------------
template<typename T>
void addMatrixMatrix(const LargeMatrix<T>& matA, const LargeMatrix<T>& matB, LargeMatrix<T>& matC)
{
  trace_p->push("addMatrixMatrix");
  if((matA.nbRows != matB.nbRows) || (matA.nbCols != matB.nbCols)  ||
      (matA.nbRows != matC.nbRows) || (matA.nbCols != matC.nbCols))
    {
      error("largematrix_mismatch_dim");
    }
  else if((matA.storage_p != matB.storage_p) ||
          (matA.storage_p->accessType() != matC.storage_p->accessType()) ||
          (matA.storage_p->storageType() != matC.storage_p->storageType()))
    {
      error("largematrix_mismatch_storage");
    }
  else if(matA.storage_p != matC.storage_p)
    {
      if(matA.values_.size() != matC.values_.size())
        {
          matC.values_.resize(matA.values_.size());
        }
      matA.storage_p->addMatrixMatrix(matA.values_, matB.values_, matC.values_);

      // Decrease number of object using matrix storage of matrix C
      matC.storage_p->objectMinus();
      matC.storage_p = matA.storage_p;

      // Increase number of object using matrix storage of matrix A (and B, of course)
      matA.storage_p->objectPlus();
    }
  else
    {
      if(matA.values_.size() != matC.values_.size())
        {
          matC.values_.resize(matA.values_.size());
        }
      matA.storage_p->addMatrixMatrix(matA.values_, matB.values_, matC.values_);
    }

  if(matA.sym != matB.sym)
    {
      matC.sym = _noSymmetry;
    }
  else if(_diagonal == matA.sym)
    {
      matC.sym = matB.sym;
    }
  else if(_diagonal == matB.sym)
    {
      matC.sym = matA.sym;
    }
  else
    {
      matC.sym = matA.sym;
    }

  trace_p->pop();
}

// ---------------------------------------------------------------------
//  Add two largeMatrix whose storage need not be the same.
//   The result matrix is based on the matrix that is most "non-symmetric"
//   This function only serves for the mode "Shift and invert" of eigenSolver.
//   Any external usage should be followed with a deletion of result matrix.
//  \param[in] matA first matrix
//  \param[in] matB second matrix
//  \return pointer to result matrix.
// ---------------------------------------------------------------------
template<typename T>
LargeMatrix<T>* addMatrixMatrixSkyline(const LargeMatrix<T>& matA, const LargeMatrix<T>& matB)
{
  LargeMatrix<T>* res = 0;
  if(matA.storagep() == matB.storagep())
    {
      res = new LargeMatrix<T>(matA);
      LargeMatrix<T>* tmpB = const_cast<LargeMatrix<T>* >(&matB);
      *res += *tmpB;
      res->toSkyline();
      return (res);
    }

  bool firstOpA = false;
  bool isAllocA = false;
  bool isAllocB = false;
  LargeMatrix<T>* pB = 0, *pA = 0;
  if(_skyline != matA.storageType())
    {
      pA = new LargeMatrix<T>(matA, true);
      pA->toSkyline();
      isAllocA = true;
    }
  else pA = const_cast<LargeMatrix<T>* >(&matA);

  if(_skyline != matB.storageType())
    {
      pB = new LargeMatrix<T>(matB, true);
      pB->toSkyline();
      isAllocB = true;
    }
  else pB = const_cast<LargeMatrix<T>* >(&matB);

  if(_dual == pA->accessType())
    {
      res = pA;
      firstOpA = true;
    }
  if(_dual == pB->accessType())
    {
      res = pB;
      firstOpA = false;
    }
  if(0 == res)
    {
      res = pA;
      firstOpA = true;
    }

  if(res == &matA)
    {
      res = new LargeMatrix<T>(matA, true);
      firstOpA = true;
    }

  if(res == &matB)
    {
      res = new LargeMatrix<T>(matB, true);
      firstOpA = false;
    }

  SkylineStorage* p = 0;
  if(firstOpA) p = static_cast<SkylineStorage*>(pB->storagep());
  else p = static_cast<SkylineStorage*>(pA->storagep());
  const std::vector<number_t>& rowPointer = p->rowPointer();
  const std::vector<number_t>& colPointer = p->colPointer();

  if(firstOpA) res->storagep()->addTwoMatrix(res->values(), res->sym, rowPointer, colPointer, pB->values(), pB->sym);
  else  res->storagep()->addTwoMatrix(res->values(), res->sym, rowPointer, colPointer, pA->values(), pA->sym);
  if((pA->values().size() != res->values().size()) || (pB->values().size() != res->values().size())) res->sym = _noSymmetry;
  if(firstOpA && isAllocB) delete pB;
  if((!firstOpA) && isAllocA) delete pA;

  return (res);
}

//  Add a scaled matrix
template<typename T>
LargeMatrix<T> addMatrixMatrix(const LargeMatrix<T>& matA, const LargeMatrix<T>& matB, T s=T(1)) //! A+s*B
{
  LargeMatrix<T> matC(matB);
  if(s!=T(1)) matC*=s;
  return matC+=matA;
}

//  Add two largeMatrix operator+
template<typename T>
LargeMatrix<T> operator+(const LargeMatrix<T>& matA, const LargeMatrix<T>& matB)
{
  LargeMatrix<T> matC(matA);
  addMatrixMatrix(matA, matB, matC);
  return matC;
}
// Specialization of operator+ of two matrices
LargeMatrix<complex_t> operator+(const LargeMatrix<real_t>& matA, const LargeMatrix<complex_t>& matB);
LargeMatrix<complex_t> operator+(const LargeMatrix<complex_t>& matA, const LargeMatrix<real_t>& matB);

//  Add two largeMatrix operator-
template<typename T>
LargeMatrix<T> operator-(const LargeMatrix<T>& matA, const LargeMatrix<T>& matB)
{
  return addMatrixMatrix(matA, matB, T(-1));
}
// Specialization of operator- of two matrices
LargeMatrix<complex_t> operator-(const LargeMatrix<real_t>& matA, const LargeMatrix<complex_t>& matB);
LargeMatrix<complex_t> operator-(const LargeMatrix<complex_t>& matA, const LargeMatrix<real_t>& matB);

/*!
   Multiple a largeMatrix with a scalar
   The result matrix will point to the same storage
   \param mat matrix
   \param v scalar
   \return result matrix which shares the same storage.
 */
template<typename T>
LargeMatrix<T> multMatrixScalar(const LargeMatrix<T>& mat, const T v)
{
  trace_p->push("multMatrixScalar");
  LargeMatrix<T> result(mat);
  typename std::vector<T>::iterator itb = result.values_.begin(), it;
  for(it = itb; it != result.values_.end(); it++)
    {
      *it *= v ;
    }
  trace_p->pop();
  return result;
}


//!  Multiple a largeMatrix with a scalar
template<typename T>
LargeMatrix<T> operator*(const LargeMatrix<T>& mat, const T v)
{
  return multMatrixScalar(mat, v);
}

//!  Multiple a largeMatrix with a scalar
template<typename T>
LargeMatrix<T> operator*(const T v, const LargeMatrix<T>& mat)
{
  return multMatrixScalar(mat, v);
}

// Specialization of multiplication largeMatrix and scalar
LargeMatrix<complex_t> operator*(const LargeMatrix<complex_t>& mat, const real_t v);
LargeMatrix<complex_t> operator*(const real_t v, const LargeMatrix<complex_t>& mat);
LargeMatrix<complex_t> operator*(const LargeMatrix<real_t>& mat, const complex_t v);
LargeMatrix<complex_t> operator*(const complex_t v, const LargeMatrix<real_t>& mat);


// ---------------------------------------------------------------------
/*! extern template product matrix*matrix
  the storage of the resulting matrix is ALWAYS dense
  there is no chance to get a sparse matrix except in case of product with
  a diagonal matrix where the storage is unchanged */
// ---------------------------------------------------------------------
template <typename SA, typename SB, typename SR>
void multMatrixMatrix(const LargeMatrix<SA>& mA, const LargeMatrix<SB>& mB, LargeMatrix<SR>& mR)
{
  if(mA.nbCols!=mB.nbRows || mA.nbColsSub!=mB.nbRowsSub)
    {
      where("multMatrixMatrix(LargeMatrix,LargeMatrix)");
      error("largematrix_mismatch_dim");
    }
  //initialisation of mR attributes
  mR.valueType_=_real;
  if(mA.valueType_==_complex || mB.valueType_==_complex) mR.valueType_=_complex;
  mR.strucType=_scalar;
  mR.nbRowsSub=mA.nbRowsSub;
  mR.nbColsSub=mB.nbColsSub;
  if(mR.nbRowsSub>1 || mR.nbColsSub>1) mR.strucType=_matrix;
  mR.nbRows=mA.nbRows;
  mR.nbCols=mB.nbCols;
  mR.sym=_noSymmetry;
  //create storage and allocate values
  if(mR.storage_p!=0) delete mR.storage_p;   //delete current storage of matrix R
  mR.storage_p=new RowDenseStorage(mR.nbRows,mR.nbCols);     //create row dense storage
  mR.storage_p->objectPlus();

  number_t s = mR.nbRows*mR.nbCols+1;
  if(Trace::traceMemory)
    {
      thePrintStream<<"LargeMatrix::multMatrixMatrix re-allocates a large matrix : "<<&mR.values_<<", "
                    <<s<<" non zeros coefficients "<<mR.dimValues();
      if(mR.storage_p!=0) thePrintStream<<", storage "<<mR.storage_p->name();
      thePrintStream<<eol<<std::flush;
    }
  mR.values_.resize(s);
  //do product
  mA.storage_p->multMatrixMatrix(mA.values_,*mB.storage_p, mB.values_,mR.values_, mA.sym, mB.sym);
}

template <typename T>
LargeMatrix<T> operator*(const LargeMatrix<T>& mA, const LargeMatrix<T>& mB)
{
  LargeMatrix<T> mR;
  multMatrixMatrix(mA,mB,mR);
  return mR;
}
//specialization
LargeMatrix<complex_t> operator*(const LargeMatrix<complex_t>&, const LargeMatrix<real_t>&);
LargeMatrix<complex_t> operator*(const LargeMatrix<real_t>&, const LargeMatrix<complex_t>&);
LargeMatrix<Matrix<complex_t> > operator*(const LargeMatrix<Matrix<complex_t> >&, const LargeMatrix<Matrix<real_t> >&);
LargeMatrix<Matrix<complex_t> > operator*(const LargeMatrix<Matrix<real_t> >&, const LargeMatrix<Matrix<complex_t> >&);

//extern template product matrix * diag matrix, diag matrix given by a vector
//when vector is a vector of d-vectors, diag matrix is assumed to be a diag matrix of d x d diag matrices
template <typename SA, typename SB, typename SR>
void multMatrixMatrix(const LargeMatrix<SA>& mA, const std::vector<SB>& mB, LargeMatrix<SR>& mR)
{
  if(mB.size()==0)
    {
      where("multMatrixMatrix(LargeMatrix,Vector)");
      error("is_void", "vector");
    }
  dimPair ds = dimsOf(mB[0]);
  if(mA.nbCols!=mB.size() || mA.nbColsSub!=ds.first)
    {
      where("multMatrixMatrix(LargeMatrix,Vector)");
      error("largematrix_mismatch_dim");
    }
  //initialisation of mR attributes
  structPair sp = Value::typeOf(mB[0]);
  mR.valueType_=_real;
  if(mA.valueType_==_complex || sp.first==_complex) mR.valueType_=_complex;
  mR.strucType=_scalar;
  mR.nbRowsSub=mA.nbRowsSub;
  mR.nbColsSub=mA.nbColsSub;
  if(mR.nbRowsSub>1 || mR.nbColsSub>1) mR.strucType=_matrix;
  mR.nbRows=mA.nbRows;
  mR.nbCols=mA.nbCols;
  mR.sym=_noSymmetry;
  //create storage and allocate values
  if(mR.storage_p!=0) delete mR.storage_p;   //delete current storage of matrix R
  mR.storage_p=mA.storage_p;                 //same storage as A
  if(Trace::traceMemory)
    {
      thePrintStream<<"LargeMatrix::multMatrixMatrix re-allocates a large matrix : "<<&mR.values_<<", "
                    <<mR.storage_p->size()<<" non zeros coefficients "<<mR.dimValues();
      if(mR.storage_p!=0) thePrintStream<<", storage "<<mR.storage_p->name();
      thePrintStream<<eol<<std::flush;
    }
  mR.values_.resize(mR.storage_p->size());
  //do product
  mA.storage_p->multMatrixDiagMatrix(mA.values_, mB,mR.values_);
}

//no operator * :  matrix * diag matrix (ambiguity with matrix * vector) , use multMarixMatrix

//extern template product diag matrix * matrix, diag matrix given by a vector
//when vector is a vector of d-vectors, diag matrix is assumed to be a diag matrix of d x d diag matrices
template <typename SA, typename SB, typename SR>
void multMatrixMatrix(const std::vector<SA>& mA, const LargeMatrix<SB>& mB, LargeMatrix<SR> & mR)
{
  if(mB.size()==0)
    {
      where("multMatrixMatrix(Vector,LargeMatrix)");
      error("is_void", "vector");
    }
  dimPair ds = dimsOf(mA[0]);
  if(mB.nbCols!=mA.size() || mB.nbColsSub!=ds.first)
    {
      where("multMatrixMatrix(Vector,LargeMatrix)");
      error("largematrix_mismatch_dim");
    }
  //initialisation of mR attributes
  structPair sp = Value::typeOf(mA[0]);
  mR.valueType_=_real;
  if(mB.valueType_==_complex || sp.first==_complex) mR.valueType_=_complex;
  mR.strucType=_scalar;
  mR.nbRowsSub=mB.nbRowsSub;
  mR.nbColsSub=mB.nbColsSub;
  if(mR.nbRowsSub>1 || mR.nbColsSub>1) mR.strucType=_matrix;
  mR.nbRows=mB.nbRows;
  mR.nbCols=mB.nbCols;
  mR.sym=_noSymmetry;
  //create storage and allocate values
  if(mR.storage_p!=0) delete mR.storage_p;   //delete current storage of matrix R
  mR.storage_p=mB.storage_p;                 //same storage as B
  if(Trace::traceMemory)
    {
      thePrintStream<<"LargeMatrix::multMatrixMatrix re-allocates a large matrix : "<<&mR.values_<<", "
                    <<mR.storage_p->size()<<" non zeros coefficients "<<mR.dimValues();
      if(mR.storage_p!=0) thePrintStream<<", storage "<<mR.storage_p->name();
      thePrintStream<<eol<<std::flush;
    }
  mR.values_.resize(mR.storage_p->size());
  //do product
  mB.storage_p->multDiagMatrixMatrix(mA,mB.values_,mR.values_);
}

//no operator * :  diag matrix * matrix (ambiguity with vector * matrix), use multMarixMatrix


// special hermitian product for vectors of form (row index, val), the second is given as a map row index -> val for better efficiency
// function used by QR algorithm
template<typename T, typename K>
T hermitianProduct(const std::vector<std::pair<number_t,T> >& u, const std::map<number_t,K>& v)
{
  typename std::vector<std::pair<number_t,T> >::const_iterator itu=u.begin();
  typename std::map<number_t,K>::const_iterator itve=v.end(), itv;
  T hp=0.;
  for(; itu!=u.end(); itu++)
    {
      itv=v.find(itu->first);
      if(itv!=itve) hp+=itu->second * conj(itv->second);
    }
  return hp;
}

// special dot product for vectors of form (row index, val), the second is given as a map row index -> val for better efficiency
// function used by QR algorithm
template<typename T, typename K>
T dotProduct(const std::vector<std::pair<number_t,T> >& u, const std::map<number_t,K>& v)
{
  typename std::vector<std::pair<number_t,T> >::const_iterator itu=u.begin();
  typename std::map<number_t,K>::const_iterator itve=v.end(), itv;
  T hp=0.;
  for(; itu!=u.end(); itu++)
    {
      itv=v.find(itu->first);
      if(itv!=itve) hp+=itu->second * itv->second;
    }
  return hp;
}

// do special linear combination u + av, used in QR algorithm
// use an (index, value) representation and eliminate near 0 components
template<typename T, typename K>
void combine(std::vector<std::pair<number_t,K> >& u, const std::map<number_t,T>& v, K a)
{
  std::vector<std::pair<number_t,K> > uv;
  uv.reserve(u.size()+v.size());
  std::set<number_t> markv;
  typename std::vector<std::pair<number_t,K> >::iterator itu=u.begin();
  typename std::map<number_t,T>::const_iterator itve=v.end(), itv;
  for(; itu!=u.end(); itu++)
    {
      itv=v.find(itu->first);
      K r= itu->second;
      if(itv!=itve)
        {
          r += a * itv->second;
          markv.insert(itu->first);
        }
      if(norm2(r)>10*theEpsilon) uv.push_back(std::pair<number_t, K>(itu->first,r));
    }
  if(markv.size()<v.size())  //add to u unmarked v components
    {
      for(itv=v.begin(); itv!=v.end(); itv++)
        if(markv.find(itv->first) == markv.end() && xlifepp::norm2(itv->second)> 10*theEpsilon)  //v component not found
          uv.push_back(std::pair<number_t,K>(itv->first, a*itv->second));
    }
  u=uv;
}

//create transposed matrix
template <typename T>
LargeMatrix<T> transpose(const LargeMatrix<T>& L)
{
  LargeMatrix<T> Lt;
  Lt.valueType_=L.valueType_;
  Lt.strucType=L.strucType;
  Lt.nbRows=L.nbCols;
  Lt.nbCols=L.nbRows;
  Lt.nbRowsSub=L.nbColsSub;
  Lt.nbColsSub=L.nbRowsSub;
  Lt.sym=L.sym;
  Lt.name=L.name+"^t";
  Lt.factorization_=_noFactorization;  //forgot factorisation properties
#ifdef XLIFEPP_WITH_UMFPACK
  Lt.umfLU_p=0;
  if(L.umfLU_p!=0)
    warning("free_warning","LargeMatrix transpose does not keep the umfpackLU object, redo factorisation if required");
#endif
  Lt.storagep()=0;
  if(L.storagep()==0)
  {
    where("transpose(LargeMatrix)");
    error("matrix_nostorage");
  }
  Lt.storagep() = L.storagep()->transpose(L.values(),Lt.values());  //build the transposed storage and fill values
  return Lt;
}

//-----------------------------------------------------------------------------------------------------------
/*!
   QR factorization of rectangular matrix pxn using Housolder method
   generic algorithm using a column algorithm with dynamic storage creation of matrix results

   mat      : matrix to factorize
   matR     : upper triangular matrix p x n stored as ColCsStorage (pointer)
   matQ     : unitary matrix of size p x p  stored as RowCsStorage (pointer)
   computeQ : if true the matrix Q is computed else not
   rhs      : right hand side vector pointer (modified), if 0 no rhs
   withColPermutation : if true, column may be permuted (more stable factorization).
            if false and algorithm fails, withColPermutation is forced to true and algorithm is rerun
   numcol_p : renumbering of columns if they have been permuted (pointer)
   stop     : iteration number when the algorithm stops: stop = nbr when algorithm exits normally

   Note 1 : to avoid unsightly rounding effect in reduction, values < 10*theEspsilon (~10^-15) are rounded to 0
   Note 2 : when there are constraint redanduncies, it may occur some average effects.
      For instance, two Dirichlet conditions at the same point (in Lagrange interpolation)
      with different right hand side produce a Dirichlet condition with a mean of right hand sides
      It may be inconvenient !
 */
//-----------------------------------------------------------------------------------------------------------

template<typename T, typename K>
void QR(const LargeMatrix<T>& mat, LargeMatrix<T>*& matR, bool computeQ, LargeMatrix<T>*& matQ, std::vector<K>* rhs,
        bool& withColPermutation, std::vector<number_t>*& numcol_p, number_t& stop)
{
  number_t nbr=mat.nbRows, nbc=mat.nbCols;
  //if(nbr > nbc) error("free_error","in ColCsStorage::QR, algorithm can no be applied when nbRows > nbCol");

  trace_p->push("QR(LargeMatrix, ...");

  //initialization of numcol
  numcol_p= new std::vector<number_t>(nbc);
  std::vector<number_t>& numcol= *numcol_p;
  std::vector<number_t>::iterator itn =numcol.begin();
  for(number_t j=0; j<nbc; j++, itn++) *itn = j;

  //construct special column representation of R matrix, for each column -> vector of pair of row index and value
  std::vector<std::vector<std::pair<number_t,T> > > cr(nbc);
  typename std::vector<std::vector<std::pair<number_t,T> > >::iterator itcr=cr.begin();
  typename std::vector<std::pair<number_t,T> >::iterator itcv;
  typename std::vector<std::pair<number_t,K> >::iterator itsv;

  number_t c=1;
  for(; itcr!=cr.end(); itcr++, c++)
    {
      std::vector<std::pair<number_t, number_t> > rowadrs=mat.getCol(c);
      std::vector<std::pair<number_t,T> > colval(rowadrs.size());
      itcv=colval.begin();
      std::vector<std::pair<number_t, number_t> >::iterator itra=rowadrs.begin();
      for(; itra!=rowadrs.end(); itra++, itcv++) *itcv= std::pair<number_t,T>(itra->first-1, mat.values()[itra->second]);
      *itcr=colval;
    }

  //construct special row representation of Q matrix, for each row -> vector of pair of col index and value
  std::vector<std::vector<std::pair<number_t,T> > > qr;
  if(computeQ)
    {
      qr.resize(nbr);
      itcr=qr.begin();
      number_t i=0;
      for(; itcr!=qr.end(); itcr++, i++) *itcr =std::vector<std::pair<number_t,T> >(1,std::pair<number_t,T>(i,T(1)));
    }

  //construct special representation of rhs vector, row index and value
  std::vector<std::pair<number_t,K> > rhsr;
  if(rhs!=0)
    {
      rhsr.resize(nbr);
      number_t i=0;
      typename std::vector<K>::iterator itv=rhs->begin();
      for(itsv=rhsr.begin(); itsv!=rhsr.end(); itsv++, itv++, i++) *itsv= std::pair<number_t,K>(i,*itv);
    }

  //to store diagonal of R matrix
  std::vector<T> diag(nbr, T(0));
  typename std::vector<T>::iterator itdiag=diag.begin();

  bool cont=true;
  number_t k=0, m=std::min(nbr,nbc);
  real_t releps=.001;

  //main loop of Housolder reduction
  //--------------------------------
  for(; k<m && cont; k++, itdiag++)
    {
      real_t nk=0;
      number_t numk=numcol[k];
      std::vector<std::pair<number_t,T> >& crk=cr[numk];
      for(itcv=crk.begin(); itcv!=crk.end(); itcv++)
        if(itcv->first>=k)
          {
            real_t nit=xlifepp::norm2(itcv->second);
            nk+=nit*nit;
          }
      real_t vmax=std::sqrt(nk);
      //locate column of largest amplitude if permutation allowed
      if(withColPermutation)
        {
          number_t jmax = k;
          for(number_t j=k+1; j< nbc; j++)  //travel column
            {
              real_t nj=0;
              std::vector<std::pair<number_t,T> >& crj=cr[numcol[j]];
              for(itcv=crj.begin(); itcv!=crj.end(); itcv++)
                if(itcv->first>=k)
                  {
                    real_t nit=xlifepp::norm2(itcv->second);
                    nj+=nit*nit;
                  }
              nj=std::sqrt(nj);
              if(nj > vmax*(1+releps))
                {
                  jmax=j;  //add criterium nj>(1+releps)vmax to avoid useless permutation
                  vmax=nj;
                }
            }
          if(jmax!=k)
            {
              numcol[k]=numcol[jmax];  //virtual permutation of columns k and j
              numcol[jmax]=numk;
            }
        }
      if(vmax< 10*theEpsilon)  cont=false;
      else
        {
          //construct Housolder vector vk
          std::map<number_t,T> vk;
          numk = numcol[k];
          for(itcv=cr[numk].begin(); itcv!=cr[numk].end(); itcv++)
            if(itcv->first >= k) vk[itcv->first]=itcv->second;  //sorted vector vk
          typename std::map<number_t,T>::iterator itvk;
          itvk=vk.find(k);
          if(itvk==vk.end()) vk[k]=-vmax;
          else itvk->second-=vmax;
          real_t nmoins=0;
          for(itvk=vk.begin(); itvk!=vk.end(); itvk++)
            {
              real_t nit=xlifepp::norm2(itvk->second);
              nmoins+=nit*nit;
            }
          vk[k]+=2*vmax;
          real_t nplus=0;
          for(itvk=vk.begin(); itvk!=vk.end(); itvk++)
            {
              real_t nit=xlifepp::norm2(itvk->second);
              nplus+=nit*nit;
            }
          if(nplus >= nmoins)
            {
              if(nmoins>10*theEpsilon) vk[k]-=2*vmax;
              else nmoins=nplus;
            }
          else
            {
              if(nplus>10*theEpsilon)  nmoins=nplus;
              else vk[k]-=2*vmax;
            };
          nmoins=std::sqrt(nmoins);
          for(itvk=vk.begin(); itvk!=vk.end(); itvk++) itvk->second/=nmoins; //Housolder vector
          std::map<number_t,T> cvk=vk;
          for(itvk=cvk.begin(); itvk!=cvk.end(); itvk++) itvk->second=conj(itvk->second); //conjugated Housolder vector

          //update R and Q from H=I-2vv*
          //itn =numcol.begin()+k;
          itn =numcol.begin();
#ifdef XLIFEPP_WITH_OMP
          #pragma omp parallel for schedule(dynamic) firstprivate(vk)
#endif
          for(number_t j=k; j<nbc; j++)
            {
              std::vector<std::pair<number_t,T> >& crj=cr[*(itn+j)];
              T psj= hermitianProduct(crj,vk);                                   //special hermitian product
              if(xlifepp::norm2(psj) > 10*theEpsilon) combine(crj, vk, -2*psj);  // do linear combination cr[rj]-2psj*vk = Hk*cr[rj]
              itcv=crj.begin();
              if(j==k)
                {
                  bool cont = true;
                  while(cont && itcv!=crj.end())
                    {
                      if(itcv->first == k)
                        {
                          *itdiag = itcv->second;
                          cont=false;
                        }
                      ++itcv;
                    }
                }
            }
          if(computeQ)
            {
              itcr=qr.begin();
              for(; itcr!=qr.end(); itcr++)
                {
                  T psi= dotProduct(*itcr,vk);                               //special dot product
                  if(norm2(psi)> 10*theEpsilon) combine(*itcr, cvk, -2*psi);// do linear combination qr[i]-2psi*conj(vk) = qr[i]*Hk^*
                }
            }
          if(rhs!=0) //update rhsr
            {
              K ps= hermitianProduct(rhsr, vk);                                //special hermitian product
              if(xlifepp::norm2(ps)> 10*theEpsilon) combine(rhsr, vk, -2*ps);  //do linear combination rhsr-2ps*vk = Hk*rhsr

            }
        } //end else
    } //end for(number_t k=0; k<m; k++)

  //end of algorithm (retry if fails)
  if(!cont && !withColPermutation) //retry QR algorithm forcing permutation of columns
    {
      warning("free_warning","QR algorithm without column permutation failed, retry forcing column permutation");
      withColPermutation=true;
      qr.clear();
      cr.clear();
      rhsr.clear();
      delete numcol_p;
      QR(mat, matR, computeQ, matQ, rhs, withColPermutation, numcol_p, stop);
    }
  else // generate final matrix
    {
      stop=m;
      if(!cont)
        {
          stop=k-1;
          warning("free_warning","QR algorithm with column permutation stopped at iteration "+tostring(stop));
        }
      //create LargeMatrix R, (we normalize R,  dividing row r by the diagonal coefficient)
      std::vector< std::vector<number_t> > rowindex(nbc);
      std::vector< std::vector<number_t> >::iterator itri=rowindex.begin();
      itn=numcol.begin();
      for(itcr=cr.begin(); itcr!=cr.end(); itcr++, itri++, itn++) //rearrange row indices
        {
          std::vector<std::pair<number_t,T> >& crc=cr[*itn];
          std::vector<number_t> ric(crc.size());
          std::vector<number_t>::iterator itk = ric.begin();
          for(itcv=crc.begin(); itcv!=crc.end(); itcv++, itk++) *itk=itcv->first +1;
          *itri = ric;
        }
      ColCsStorage* css=new ColCsStorage(nbr, nbc, rowindex);  //create col storage
      matR=new LargeMatrix<T>(mat.valueType_, mat.strucType, nbr, nbc, _noSymmetry, mat.nbRowsSub, mat.nbColsSub, T(0), "R", css);
      typename std::vector<T>::iterator itval=matR->values().begin()+1;
      itn=numcol.begin();
      for(itcr=cr.begin(); itcr!=cr.end(); itcr++, itn++)  //sequential copy of values
        {
          std::vector<std::pair<number_t,T> >& crc=cr[*itn];
          for(itcv=crc.begin(); itcv!=crc.end(); itcv++, itval++)
            {
              *itval = itcv->second;
              number_t r=itcv->first;
              T d=diag[r];
              if(d!=T(0)) *itval /=d;
            }
        }
      //create LargeMatrix Q if required
      matQ=0;
      if(computeQ)
        {
          std::vector< std::vector<number_t> > colindex(nbr);
          std::vector< std::vector<number_t> >::iterator itci=colindex.begin();
          for(itcr=qr.begin(); itcr!=qr.end(); itcr++, itci++) //rearrange row indices
            {
              std::vector<std::pair<number_t,T> >& qri=*itcr;
              std::vector<number_t> ric(qri.size());
              itn=ric.begin();
              for(itcv=qri.begin(); itcv!=qri.end(); itcv++, itn++) *itn=itcv->first +1;
              *itci = ric;
            }
          RowCsStorage* rss=new RowCsStorage(nbr, nbr, colindex);  //create row storage
          matQ=new LargeMatrix<T>(mat.valueType_, mat.strucType, nbr, nbr, _noSymmetry, mat.nbRowsSub, mat.nbColsSub, T(0), "Q", rss);
          itval=matQ->values().begin()+1;
          for(itcr=qr.begin(); itcr!=qr.end(); itcr++)  //sequential copy of values
            {
              std::vector<std::pair<number_t,T> >& qri=*itcr;
              for(itcv=qri.begin(); itcv!=qri.end(); itcv++, itval++) *itval = itcv->second;
            }
        }
      if(rhs!=0) //update rhs
        {
          rhs->clear();
          rhs->resize(nbr, K(0));
          typename std::vector<K>::iterator itv=rhs->begin();
          for(itsv=rhsr.begin(); itsv!=rhsr.end(); itsv++)
            {
              number_t r=itsv->first;
              T d=diag[r];
              if(d!=T(0)) * (itv + r) = itsv->second / d;  // normalization
              else * (itv + r) = itsv->second;
            }
        }

      //finalization
      if(!withColPermutation)
        {
          delete numcol_p;
          numcol_p=0;
        }
    }

  trace_p->pop();
}
/*!
   solve upper triangular system pxp with diag = Id with a list of right hand sides
   each right hand side is a vector of pairs (row index, value)
   using a column algorithm with dynamic storage update
          p
      | 1 x ... x |
      | 0 1 ... x |
 at = |     ...   |
      | 0   0 1 x |
      | 0     0 1 |

   mat      : upper triangular matrix
   rhss     : list of right hand sides

   rhss is overwritten !

*/
template<typename T, typename K>
void QRSolve(const LargeMatrix<T>& mat, std::vector<std::vector<std::pair<number_t,K> > >& rhss)
{
  number_t p=mat.nbRows, k = p;
  typename std::vector<std::vector<std::pair<number_t,K> > >::iterator itr;
  typename std::vector<std::pair<number_t,K> >::iterator itp;
  const std::vector<T>& val=mat.values();
  for(; k>1; k--)
    {
      std::vector<std::pair<number_t, number_t> > colpadr=mat.getCol(k, 1, k-1);
      if(colpadr.size()>0)
        {
          std::vector<std::pair<number_t, number_t> >::iterator ita=colpadr.begin();
          std::map<number_t,T> mapval;
          for(; ita!=colpadr.end(); ita++) mapval.insert(std::make_pair(ita->first-1,val[ita->second]));
          for(itr=rhss.begin(); itr!=rhss.end(); itr++)
            {
              itp=itr->begin();
              bool cont=true;
              while(itp!=itr->end() && cont)   //to be improved (sequential search)
                {
                  if(itp->first == k-1)
                    {
                      combine(*itr, mapval,-itp->second);
                      cont=false;
                    }
                  itp++;
                }
            }
        }
    }
}

// solve upper triangular system mat * x = b with some matrix columns and one vector as right hand sides
// diagonal coefficient of mat are assumed to be 1
// create a general right hand sides with storage well adapted to dynamic column size computation used in previous function
// Note : this particular form has been developed in the context of management of essential conditions where matrices are very sparse
// WARNING : this algorithm does not work yet if template T and K are not the same !!! TODO in future
template<typename T, typename K>
void QRSolve(const LargeMatrix<T>& mat, LargeMatrix<T>* mat2, std::vector<K>* rhs)
{
  if(mat2==0 && rhs==0) return;
  trace_p->push("QRSolve(LargeMatrix, ...");

  number_t nr=0, nbr=mat.nbRows, nbc=0;
  if(rhs!=0) nr+=1;
  if(mat2!=0)
    {
      nbc=mat2->nbCols;
      nr+=nbc;
    }
  std::vector<std::vector<std::pair<number_t,K> > > rhss;
  rhss.resize(nr);
  typename std::vector<std::vector<std::pair<number_t,K> > >::iterator itr=rhss.begin();
  typename std::vector<std::pair<number_t, K> >::iterator its;
  typename std::vector<std::pair<number_t, T> >::iterator itv;

  if(mat2!=0)
    {
      std::vector<T>& val=mat2->values();
      for(number_t c=1; c<=nbc; c++, itr++)
        {
          std::vector<std::pair<number_t, number_t> > colpadr=mat2->getCol(c);
          std::vector<std::pair<number_t, number_t> >::iterator ita=colpadr.begin();
          itr->resize(colpadr.size());
          itv=itr->begin();
          for(; ita!=colpadr.end(); ita++, itv++) *itv = std::make_pair(ita->first-1,val[ita->second]);
        }
    }

  if(rhs!=0)
    {
      itr->resize(nbr);
      itv=itr->begin();
      typename std::vector<T>::iterator its=rhs->begin();
      for(number_t r=0; r<nbr; r++, itv++, its++) *itv = std::make_pair(r,*its);
    }

  QRSolve(mat,rhss);   //call general QRSolver

  if(mat2!=0)   //reconstruct mat2
    {
      mat2->clear();
      std::vector< std::vector<number_t> > rowindex(nbc);
      std::vector< std::vector<number_t> >::iterator itri=rowindex.begin();
      itr=rhss.begin();
      for(number_t c=0; c<nbc; c++, itr++, itri++)
        {
          itri->resize(itr->size());
          std::vector<number_t>::iterator itk = itri->begin();
          for(itv=itr->begin(); itv!=itr->end(); itv++, itk++) *itk=itv->first +1;
        }
      ColCsStorage* css=new ColCsStorage(nbr, nbc, rowindex);  //create col storage
      mat2->init(css, T(0), _noSymmetry);
      typename std::vector<T>::iterator itval=mat2->values().begin()+1;
      itr=rhss.begin();
      for(number_t c=0; c<nbc; c++, itr++)
        for(itv=itr->begin(); itv!=itr->end(); itv++, itval++) *itval=itv->second;
    }

  if(rhs!=0) //reconstruct rhs
    {
      rhs->clear();
      rhs->resize(nbr,T(0));
      for(its=itr->begin(); its!=itr->end(); its++) *(rhs->begin()+ its->first) =its->second;
    }
  trace_p->pop();
}

//Gauss solver using algorithms defined in storage classes
template<typename T>
void gaussSolve(LargeMatrix<T>& mat, std::vector<std::vector<T> >& rhss)
{
  mat.storagep()->gaussSolver(mat.values(),rhss);
}

template<typename T>
void gaussSolve(LargeMatrix<T>& mat, std::vector<T>& rhs)
{
  mat.storagep()->gaussSolver(mat.values(),rhs);
}

//-----------------------------------------------------------------------------------------------------------------------------
/*!
  Resolve an eigen problem with block davidson method
  \param[in] pA Pointer to large matrix A of the eigen problem A*X = lambda*B*X. It MUSTN'T be NULL
  \param[in] pB Pointer to large matrix B of the eigen problem A*X = lambda*B*X. If it's NULL. Eigen problem is standard: A*X=lambda*X
  \param[in,out] res result
  \param[in] nev The number of eigen values and eigenSolver searched for
  \param[in] tol Tolerance
  \param[in] which Specification of which eigen values are returned. Largest-LM, Smallest-SM,
  \param[in] isInverted true if inverted matrix
  \param[in] fac factorization type
  \param[in] isShift true if eigen problem with shift or not
*/
template<typename ST>
number_t eigenDavidsonSolve(const LargeMatrix<ST>* pA, const LargeMatrix<ST>* pB, std::vector<std::pair<complex_t,Vector<complex_t> > >& res,
                            number_t nev, real_t tol, string_t which, bool isInverted, FactorizationType fac, bool isShift)
{
  if(0 == pA)
    {
      where("eigenDavidsonSolve(...)");
      error("is_void", "matrix A");
    }
  if((0 != pB) && ((_symmetric != pB->sym) && (_selfAdjoint != pB->sym)))
    {
      where("eigenDavidsonSolve(...)");
      error("largematrix_mismatch_sym");
    }
  if((0 != pA) && (0 != pB) && (pA->nbCols != pB->nbCols))
    {
      where("eigenDavidsonSolve(...)");
      error("largematrix_mismatch_dim");
    }
  if(nev > pA->nbCols)
    {
      where("eigenDavidsonSolve(...)");
      error("is_greater", nev, pA->nbCols);
    }

  typedef NumTraits<ST> SCT;
  typedef MultiVec<ST>  MV;
  typedef Operator<ST>  OP;

  string_t ortho("SVQB");
  bool locking = true;
  bool insitu = false;
  MsgEigenType verbosity = _errorsEigen;

  number_t probSize = pA->nbCols;
  number_t maxNumBlocks = 5;
  int_t maxLocks = nev;
  int_t maxRestarts = 200;

  // Compute automatically nev and blockSize as well as numblock
  number_t blockSize = nev;
  number_t numBlocks = 2;
  while((blockSize*numBlocks) > probSize) --blockSize;
  for(number_t i = 3 ; i <= maxNumBlocks; ++i)
    {
      if((blockSize*i) > probSize) break;
      numBlocks = i;
    }

  while(((blockSize+maxLocks) > nev) && ((blockSize*numBlocks+maxLocks) > probSize)) --maxLocks;
  if(maxLocks < 0)
    {
      maxLocks = 0;
      locking = false;
    }
  while((numBlocks>2) && ((blockSize*numBlocks+maxLocks) > probSize)) --numBlocks;
  if((blockSize*numBlocks+maxLocks) > probSize)
    {
      nev = static_cast<int>(std::floor(probSize/numBlocks));
      blockSize = nev;
      maxLocks  = probSize - (blockSize*numBlocks);
      string_t msg = "ATTENTION! LargeMatrix::eigenDavidsonSolver. Because the matrix is too small. Only ";
      msg += tostring(nev) + " eigenvalues can be retrieved! \n";
      warning("free_warning",msg);
    }

  // Create parameter list to pass into the solver manager
  Parameters pl;
  pl.add("Verbosity", (int_t)verbosity);
  pl.add("Which", which);
  pl.add("Orthogonalization", ortho);
  pl.add("Block Size", blockSize);
  pl.add("Num Blocks", numBlocks);
  pl.add("Maximum Restarts", maxRestarts);
  pl.add("Convergence Tolerance", tol);
  pl.add("Use Locking", locking);
  pl.add("Locking Tolerance", tol/10);
  pl.add("In Situ Restarting", insitu);
  pl.add("Max Locked", maxLocks);

  SmartPtr<const OP> A = _smPtrNull, B = _smPtrNull;
  if(0 != pB)
    {
      if(isShift)
        {
          pl.add("Which", string_t("LM"));
          A = new LargeMatrixAdapterInverseGeneralized<LargeMatrix<ST>,ST>(pA, pB, fac);
        }
      else
        {
          if("SM" == which)
            {
              pl.add("Which", string_t("LM"));
              A = new LargeMatrixAdapterInverseGeneralized<LargeMatrix<ST>,ST>(pA, pB, fac);
            }
          else
            {
              A = new LargeMatrixAdapter<LargeMatrix<ST>,ST>(pA);
              B = new LargeMatrixAdapter<LargeMatrix<ST>,ST>(pB);
            }
        }
    }
  else
    {
      if(isShift)
        {
          A = new LargeMatrixAdapterInverse<LargeMatrix<ST>,ST>(pA, fac);
        }
      else
        {
          if("SM" == which)
            {
              pl.add("Which", string_t("LM"));
              A = new LargeMatrixAdapterInverse<LargeMatrix<ST>,ST>(pA, fac);
            }
          else
            {
              A = new LargeMatrixAdapter<LargeMatrix<ST>,ST>(pA);
            }
        }
    }

  SmartPtr<MV> ivec = SmartPtr<MultiVecAdapter<ST> >(new MultiVecAdapter<ST>(pA->nbCols, blockSize));
  ivec->mvRandom();

  // Define eigen problem
  SmartPtr<EigenProblem<ST,MV,OP> > problem;
  if(_smPtrNull != B) problem = SmartPtr<EigenProblem<ST,MV,OP> >(new EigenProblem<ST,MV,OP>(A,B,ivec));
  else problem = SmartPtr<EigenProblem<ST,MV,OP> >(new EigenProblem<ST,MV,OP>(A,ivec));

  // Inform the eigenproblem that the operator cK is symmetric
  problem->setHermitian(true);

  // Set the number of eigenvalues requested
  problem->setNEV(nev);

  // Inform the eigenproblem that you are done passing it information
  bool boolret = problem->setProblem();
  if(!boolret) error("eigen_eigenproblem", "");

  BlockDavidsonSolMgr<ST,MV,OP> bdSolverMan(problem, pl);
  ComputationInfo returnCode = bdSolverMan.solve();

  number_t numev = 0;
  if(_success == returnCode)
    {
      EigenSolverSolution<ST,MV> sol = problem->getSolution();
      SmartPtr<MV> evecs = sol.evecs;
      numev = sol.numVecs;
      res.resize(numev);
      for(number_t i = 0; i<numev; i++)
        {
          std::vector<complex_t> eVec(ivec->getVecLength());
          for(number_t j = 0; j < ivec->getVecLength(); ++j)
            {
              if(0 == sol.index[i]) eVec[j] = (complex_t(SCT::real((*evecs)(j,i)),0));
              else if(1 == sol.index[i]) eVec[j] = (complex_t(SCT::real((*evecs)(j,i)),SCT::real((*evecs)(j,i+1))));
              else eVec[i] = (complex_t(SCT::real((*evecs)(j,i-1)),SCT::real((*evecs)(j,i))));
            }
          res.at(i) = std::make_pair(complex_t(sol.evals[i].realpart,sol.evals[i].imagpart), eVec);
        }
    }

  return (numev);
}

//-----------------------------------------------------------------------------------------------------------------------------
/*!
  Resolve an eigen problem with block Krylov-Schur method
 \param[in] pA Pointer to large matrix A of the eigen problem A*X = lambda*B*X. It MUSTN'T be NULL
 \param[in] pB Pointer to large matrix B of the eigen problem A*X = lambda*B*X. If it's NULL. Eigen problem is standard: A*X=lambda*X
 \param[in,out] res result
 \param[in] nev The number of eigen values and eigenSolver searched for
 \param[in] tol Tolerance
 \param[in] which Specification of which eigen values are returned. Largest-LM, Smallest-SM,
 \param[in] isInverted true if inverted matrix
 \param[in] fac factorization type
 \param[in] isShift true if eigen problem with shift or not
*/
template<typename ST>
number_t eigenKrylovSchurSolve(const LargeMatrix<ST>* pA, const LargeMatrix<ST>* pB, std::vector<std::pair<complex_t,Vector<complex_t> > >& res,
                               number_t nev, real_t tol, string_t which, bool isInverted, FactorizationType fac, bool isShift)
{
  if (0 == pA)
  {
    where("eigenKrylovSchurSolve(...)");
    error("null_pointer", "matrix A");
  }
  if (nev > pA->nbCols)
  {
    where("eigenKrylovSchurSolve(...)");
    error("is_greater", nev, pA->nbCols);
  }
  typedef NumTraits<ST> SCT;
  typedef MultiVec<ST>  MV;
  typedef Operator<ST>  OP;

  string_t ortho("SVQB");
  bool insitu = false;
  MsgEigenType verbosity = _errorsEigen;

  number_t probSize = pA->nbCols;
  number_t maxNumBlocks = (nev < 3) ? 3 : nev;
  number_t maxRestarts = 250;

  //
  // Compute automatically nev and blockSize as well as numblock
  //
  number_t blockSize = 1;
  number_t numBlocks = (nev/blockSize < 3) ? 3 : (static_cast<int>(std::ceil(nev/blockSize)));
  number_t i = 0;
  while(i<(maxNumBlocks-3))
    {
      if((blockSize*numBlocks) > probSize) break;
      ++numBlocks;
      ++i;
    }
  if((blockSize*numBlocks) > probSize)
    {
      while((blockSize >1) && ((blockSize*numBlocks) > probSize)) --blockSize;
      while((numBlocks >3) && ((blockSize*numBlocks+blockSize) > probSize)) --numBlocks;
      if((probSize-1) == numBlocks) numBlocks -= blockSize;
      nev = numBlocks-1; //blockSize;
      string_t msg = "ATTENTION! LargeMatrix::eigenKrylovSchurSolver. Because the matrix is too small. Only ";
      msg += tostring(nev) + " eigenvalues can be retrieved! \n";
      warning("free_warning",msg);
    }

  // Create parameter list to pass into the solver manager
  Parameters pl;
  pl.add("Verbosity", (int_t)verbosity);
  pl.add("Which", which);
  pl.add("Orthogonalization", ortho);
  pl.add("Block Size", (int_t)blockSize);
  pl.add("Num Blocks", (int_t)numBlocks);
  pl.add("Maximum Restarts", (int_t)maxRestarts);
  pl.add("Convergence Tolerance", tol);
  pl.add("Locking Tolerance", tol/10);
  pl.add("In Situ Restarting", insitu);

//    LargeMatrix<ST>* pcB = 0;
  SmartPtr<const OP> A = _smPtrNull, B = _smPtrNull;
  if(0 == pB)
    {
      if(isShift)
        {
          pl.add("Which", string_t("LM"));
          A = new LargeMatrixAdapterInverse<LargeMatrix<ST>,ST>(pA, fac);
        }
      else
        {
          if("SM" == which)
            {
              A = new LargeMatrixAdapterInverse<LargeMatrix<ST>,ST>(pA, fac);
              pl.add("Which", string_t("LM"));
            }
          else
            {
              A = new LargeMatrixAdapter<LargeMatrix<ST>,ST>(pA);
            }
        }
    }
  else
    {
      if(isShift)
        {
          pl.add("Which", string_t("LM"));
          A = new LargeMatrixAdapterInverseGeneralized<LargeMatrix<ST>,ST>(pA, pB, fac);
          B = new LargeMatrixAdapter<LargeMatrix<ST>,ST>(pB);
        }
      else
        {
          if("SM" == which)
            {
              pl.add("Which", string_t("LM"));
              A = new LargeMatrixAdapterInverseGeneralized<LargeMatrix<ST>,ST>(pA, pB, fac);  // A^-1*B*X=lambda*X
              B = new LargeMatrixAdapter<LargeMatrix<ST>,ST>(pB);
            }
          else
            {
              A = new LargeMatrixAdapterInverseGeneralized<LargeMatrix<ST>,ST>(pB, pA, fac);  // B^-1*A*X=lambda*X
              B = new LargeMatrixAdapter<LargeMatrix<ST>,ST>(pA);
            }
        }
    }

  SmartPtr<MV> ivec = SmartPtr<MultiVecAdapter<ST> >(new MultiVecAdapter<ST>(pA->nbCols, blockSize));
  ivec->mvRandom();

  // Define eigen problem
  SmartPtr<EigenProblem<ST,MV,OP> > problem;
  if(0 != pB)  problem = SmartPtr<EigenProblem<ST,MV,OP> >(new EigenProblem<ST,MV,OP>(A,B,ivec));
  else problem = SmartPtr<EigenProblem<ST,MV,OP> >(new EigenProblem<ST,MV,OP>(A,ivec));

  // Inform the eigenproblem that the operator cK is symmetric or non-symmetric
  // By default, it's set to non-symmetric
  problem->setHermitian(false);

  // Set the number of eigenvalues requested
  problem->setNEV(nev);

  // Inform the eigenproblem that you are done passing it information
  bool boolret = problem->setProblem();
  if(!boolret) error("eigen_eigenproblem", "");

  BlockKrylovSchurSolMgr<ST,MV,OP> krSolverMan(problem, pl);
  ComputationInfo returnCode = krSolverMan.solve();

  number_t numev = 0;
  if(_success == returnCode)
    {
      EigenSolverSolution<ST,MV> sol = problem->getSolution();
      SmartPtr<MV> evecs = sol.evecs;
      numev = sol.numVecs;
      res.resize(numev);
      for(number_t i = 0; i<numev; i++)
        {
          std::vector<complex_t> eVec(ivec->getVecLength());
          for(number_t j = 0; j < ivec->getVecLength(); ++j)
            {
              if(0 == sol.index[i]) eVec[j] = (complex_t(SCT::real((*evecs)(j,i)),0));
              else if(1 == sol.index[i]) eVec[j] = (complex_t(SCT::real((*evecs)(j,i)),SCT::real((*evecs)(j,i+1))));
              else eVec[i] = (complex_t(SCT::real((*evecs)(j,i-1)),SCT::real((*evecs)(j,i))));
            }
          res.at(i) = std::make_pair(complex_t(sol.evals[i].realpart,sol.evals[i].imagpart), eVec);
        }
    }

  return (numev);
}

} // end of namespace xlifepp

#endif // LARGE_MATRIX_HPP