src/largeMatrix/MatrixEntry.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
  \file MatrixEntry.cpp
  \author E. Lunéville
  \since 01 nov 2013
  \date 01 nov 2013

  \brief Implementation of xlifepp::MatrixEntry class member functions and related utilities
*/

#include "MatrixEntry.hpp"
#include "utils.h"
#ifdef XLIFEPP_WITH_UMFPACK
#include "umfpackSupport.h"
#endif
namespace xlifepp
{

//===============================================================================
//member functions of MatrixEntry class
//===============================================================================

//constructors
MatrixEntry::MatrixEntry()
{
   rEntries_p = 0; cEntries_p = 0; rmEntries_p = 0; cmEntries_p = 0;
   valueType_ = _real; strucType_ = _scalar; nbOfComponents = dimPair(1, 1);
}

MatrixEntry::MatrixEntry(ValueType vt, StrucType st, MatrixStorage* ms_p, dimPair valdim, SymType sy)
{
   rEntries_p = 0; cEntries_p = 0; rmEntries_p = 0; cmEntries_p = 0;
   valueType_ = vt; strucType_ = st; nbOfComponents = valdim;

   if(strucType_ == _scalar && valdim != dimPair(1, 1))
   {
      error("matrixentry_incoherent_dim", words("structure", _scalar), valdim.first, valdim.second);
   }
   if(strucType_ == _matrix && valdim == dimPair(1, 1))
   {
      error("matrixentry_incoherent_dim", words("structure", _matrix), 1, 1);
   }

   switch(st)
   {
      case _scalar :
         switch(vt)
         {
            case _real :    rEntries_p = new LargeMatrix<real_t>(ms_p, sy); return;
            case _complex : cEntries_p = new LargeMatrix<complex_t>(ms_p, sy); return;
            default : error("matrixentry_abnormal_type", words("value", vt));

         }
      case _matrix :
         switch(vt)
         {
            case _real :    rmEntries_p = new LargeMatrix<Matrix<real_t> >(ms_p, Matrix<real_t>(valdim.first, valdim.second, 0.), sy); return;
            case _complex : cmEntries_p = new LargeMatrix<Matrix<complex_t> >(ms_p, Matrix<complex_t>(valdim.first, valdim.second, 0.), sy); return;
            default : error("matrixentry_abnormal_type", words("value", vt));
         }
      default : error("matrixentry_novector");
   }
}

//! construct special matrix (constant diagonal)
MatrixEntry::MatrixEntry(SpecialMatrix spm, StorageType st, AccessType at, number_t nbr, number_t nbc, const real_t& a)
{
   valueType_ = _real; strucType_ = _scalar; nbOfComponents = dimPair(1, 1);
   cEntries_p = 0; rmEntries_p = 0; cmEntries_p = 0;
   rEntries_p=new LargeMatrix<real_t>(spm, st, at , nbr, nbc, a);
}

MatrixEntry::MatrixEntry(SpecialMatrix spm, StorageType st, AccessType at, number_t nbr, number_t nbc, const complex_t& a)
{
   valueType_ = _complex; strucType_ = _scalar; nbOfComponents = dimPair(1, 1);
   rEntries_p = 0; rmEntries_p = 0; cmEntries_p = 0;
   cEntries_p=new LargeMatrix<complex_t>(spm, st, at , nbr, nbc, a);
}

MatrixEntry::MatrixEntry(SpecialMatrix spm, StorageType st, AccessType at, number_t nbr, number_t nbc, const Matrix<real_t>& a)
{
   valueType_ = _real; strucType_ = _matrix; nbOfComponents = dimPair(a.numberOfRows(), a.numberOfColumns());
   rEntries_p = 0; cEntries_p = 0; cmEntries_p = 0;
   rmEntries_p=new LargeMatrix<Matrix<real_t> >(spm, st, at , nbr, nbc, a);
}

MatrixEntry::MatrixEntry(SpecialMatrix spm, StorageType st, AccessType at, number_t nbr, number_t nbc, const Matrix<complex_t>& a)
{
   valueType_ = _complex; strucType_ = _matrix; nbOfComponents = dimPair(a.numberOfRows(), a.numberOfColumns());
   rEntries_p = 0; cEntries_p = 0; rmEntries_p = 0;
   cmEntries_p=new LargeMatrix<Matrix<complex_t> >(spm, st, at , nbr, nbc, a);
}

//copy constructor and assignment operator, it is a full copy !!!
MatrixEntry::MatrixEntry(const MatrixEntry& mat, bool storageCopy)
{
   copy(mat, storageCopy);
}

MatrixEntry& MatrixEntry::operator=(const MatrixEntry& mat)
{
   clear();
   copy(mat);   //does not copy the storage
   return *this;
}

// deallocate memory used by a MatrixEntry
void MatrixEntry::clear()
{
   if(rEntries_p != 0) delete rEntries_p;
   if(cEntries_p != 0) delete cEntries_p;
   if(rmEntries_p != 0) delete rmEntries_p;
   if(cmEntries_p != 0) delete cmEntries_p;
   rEntries_p = 0; cEntries_p = 0; rmEntries_p = 0; cmEntries_p = 0;
}

//copy tool of a MatrixEntry, object has to be clean before
void MatrixEntry::copy(const MatrixEntry& mat, bool storageCopy)
{
   valueType_ = mat.valueType_;
   strucType_ = mat.strucType_;
   nbOfComponents = mat.nbOfComponents;
   rEntries_p = 0; cEntries_p = 0; rmEntries_p = 0; cmEntries_p = 0;
   switch(strucType_)
   {
      case _scalar :
         switch(valueType_)
         {
            case _real :    if(mat.rEntries_p!=0) rEntries_p = new LargeMatrix<real_t>(*mat.rEntries_p, storageCopy); return;
            case _complex : if(mat.cEntries_p!=0) cEntries_p = new LargeMatrix<complex_t>(*mat.cEntries_p, storageCopy); return;
            default : error("matrixentry_abnormal_type", words("value", valueType_));
         }
      case _matrix :
         switch(valueType_)
         {
            case _real :    if(mat.rmEntries_p!=0) rmEntries_p = new LargeMatrix<Matrix<real_t> >(*mat.rmEntries_p, storageCopy); return;
            case _complex : if(mat.cmEntries_p!=0) cmEntries_p = new LargeMatrix<Matrix<complex_t> >(*mat.cmEntries_p, storageCopy); return;
            default : error("matrixentry_abnormal_type", words("value", valueType_));
         }
      default : error("matrixentry_novector");
   }
}

//destructor
MatrixEntry::~MatrixEntry()
{
   if(rEntries_p != 0) delete rEntries_p;
   if(cEntries_p != 0) delete cEntries_p;
   if(rmEntries_p != 0) delete rmEntries_p;
   if(cmEntries_p != 0) delete cmEntries_p;
}

// print entry
void MatrixEntry::print(std::ostream& out) const
{
   if(rEntries_p != 0)  {rEntries_p->print(out); return;}
   if(cEntries_p != 0)  {cEntries_p->print(out); return;}
   if(rmEntries_p != 0) {rmEntries_p->print(out); return;}
   if(cmEntries_p != 0) {cmEntries_p->print(out); return;}
   out<<" void matrix"<<eol;
}

// return the number of rows
number_t MatrixEntry::nbOfRows() const
{
   if(rEntries_p != 0)  return rEntries_p->nbRows;
   if(cEntries_p != 0)  return cEntries_p->nbRows;
   if(rmEntries_p != 0) return rmEntries_p->nbRows;
   if(cmEntries_p != 0) return cmEntries_p->nbRows;
   return 0;
}

// return the number of columns
number_t MatrixEntry::nbOfCols() const
{
   if(rEntries_p != 0)  return rEntries_p->nbCols;
   if(cEntries_p != 0)  return cEntries_p->nbCols;
   if(rmEntries_p != 0) return rmEntries_p->nbCols;
   if(cmEntries_p != 0) return cmEntries_p->nbCols;
   return 0;
}

// set the number of columns
void MatrixEntry::setNbOfCols(number_t n)
{
   if(rEntries_p != 0)  rEntries_p->nbCols=n;
   if(cEntries_p != 0)  cEntries_p->nbCols=n;
   if(rmEntries_p != 0) rmEntries_p->nbCols=n;
   if(cmEntries_p != 0) cmEntries_p->nbCols=n;
}

// set the number of rows
void MatrixEntry::setNbOfRows(number_t n)
{
   if(rEntries_p != 0)  rEntries_p->nbRows=n;
   if(cEntries_p != 0)  cEntries_p->nbRows=n;
   if(rmEntries_p != 0) rmEntries_p->nbRows=n;
   if(cmEntries_p != 0) cmEntries_p->nbRows=n;
}

//storage pointer, storage type and access type
StorageType MatrixEntry::storageType() const
{
   if(rEntries_p != 0)  return rEntries_p->storageType();
   if(cEntries_p != 0)  return cEntries_p->storageType();
   if(rmEntries_p != 0) return rmEntries_p->storageType();
   if(cmEntries_p != 0) return cmEntries_p->storageType();
   return _noStorage;
}

AccessType MatrixEntry::accessType() const
{
   if(rEntries_p != 0)  return rEntries_p->accessType();
   if(cEntries_p != 0)  return cEntries_p->accessType();
   if(rmEntries_p != 0) return rmEntries_p->accessType();
   if(cmEntries_p != 0) return cmEntries_p->accessType();
   return _noAccess;
}

MatrixStorage* MatrixEntry::storagep() const
{
   if(rEntries_p != 0)  return rEntries_p->storagep();
   if(cEntries_p != 0)  return cEntries_p->storagep();
   if(rmEntries_p != 0) return rmEntries_p->storagep();
   if(cmEntries_p != 0) return cmEntries_p->storagep();
   return 0;
}

MatrixStorage* MatrixEntry::storagep()
{
   if(rEntries_p != 0)  return rEntries_p->storagep();
   if(cEntries_p != 0)  return cEntries_p->storagep();
   if(rmEntries_p != 0) return rmEntries_p->storagep();
   if(cmEntries_p != 0) return cmEntries_p->storagep();
   return 0;
}


SymType MatrixEntry::symmetry() const
{
   if(rEntries_p != 0)  return rEntries_p->sym;
   if(cEntries_p != 0)  return cEntries_p->sym;
   if(rmEntries_p != 0) return rmEntries_p->sym;
   if(cmEntries_p != 0) return cmEntries_p->sym;
   return _noSymmetry;
}

SymType& MatrixEntry::symmetry()
{
  if(rEntries_p != 0)  return rEntries_p->sym;
  if(cEntries_p != 0)  return cEntries_p->sym;
  if(rmEntries_p != 0) return rmEntries_p->sym;
  if(cmEntries_p != 0) return cmEntries_p->sym;
  where("MatrixEntry::symmetry");
  error("null_pointer","xxEntries_p");
  return rEntries_p->sym; // dummy return
}
//! return a reference to LargeMatrix object
template <>
LargeMatrix<real_t>& MatrixEntry::getLargeMatrix() const
{
  if (rEntries_p!=0) return *rEntries_p;
  where("MatrixEntry::getLargeMatrix()");
  error("null_pointer","rEntries_p");
  return *rEntries_p;  //fake return
}

template <>
LargeMatrix<complex_t>& MatrixEntry::getLargeMatrix() const
{
  if (cEntries_p!=0) return *cEntries_p;
  where("MatrixEntry::getLargeMatrix()");
  error("null_pointer","cEntries_p");
  return *cEntries_p;  //fake return
}

template <>
LargeMatrix<Matrix<real_t> >& MatrixEntry::getLargeMatrix() const
{
  if(rmEntries_p!=0) return *rmEntries_p;
  where("MatrixEntry::getLargeMatrix()");
  error("null_pointer","rmEntries_p");
  return *rmEntries_p;  //fake return
}

template <>
LargeMatrix<Matrix<complex_t> >& MatrixEntry::getLargeMatrix() const
{
  if(cmEntries_p!=0) return *cmEntries_p;
  where("MatrixEntry::getLargeMatrix()");
  error("null_pointer","cmEntries_p");
  return *cmEntries_p;  //fake return
}

// return a row or a col as VectorEntry
//  rc : row or col index
//  at : row acces (_row) or col access (_col)
//  ### only available for scalar MatrixEntry ###
VectorEntry MatrixEntry::getRowCol(number_t rc, AccessType at) const
{
  //get (col/row index, adress) pairs of row/col rc
  std::vector<std::pair<number_t, number_t> > indadd;
  number_t n=nbOfRows();
  if (at==_row) n=nbOfCols();
  SymType sym=symmetry();
  if (rEntries_p!=0)
  {
    if (at==_row) indadd=rEntries_p->getRow(rc); else indadd=rEntries_p->getCol(rc);
    VectorEntry v(_real,_scalar, n);
    std::vector<real_t>::iterator itv=v.rEntries_p->begin();
    std::vector<std::pair<number_t, number_t> >::iterator itp=indadd.begin();
    std::vector<real_t>& mat = rEntries_p->values();
    switch (sym)
    {
      case _skewSymmetric:
      case _skewAdjoint:
        for (;itp!=indadd.end(); ++itp)
        {
            number_t cr = itp->first;
            if ((at==_row && rc > cr) || (at==_col && rc < cr)) *(itv+(cr-1))= -mat[itp->second];
            else *(itv+(itp->first-1))=mat[itp->second];
        }
        break;
      default :
        for (;itp!=indadd.end(); ++itp) *(itv+(itp->first-1))=mat[itp->second];
    }
    return v;
  }
  if (cEntries_p!=0)
  {
    if (at==_row) indadd=rEntries_p->getRow(rc); else indadd=rEntries_p->getCol(rc);
    VectorEntry v(_complex,_scalar,n);
    std::vector<complex_t>::iterator itv=v.cEntries_p->begin();
    std::vector<std::pair<number_t, number_t> >::iterator itp=indadd.begin();
    std::vector<complex_t>& mat = cEntries_p->values();
    switch (sym)
    {
      case _skewSymmetric:
        for (;itp!=indadd.end(); ++itp)
        {
          number_t cr = itp->first;
          if ((at==_row && rc > cr) || (at==_col && rc < cr)) *(itv+(cr-1))= -mat[itp->second];
          else *(itv+(itp->first-1))=mat[itp->second];
        }
        break;
      case _skewAdjoint:
        for (;itp!=indadd.end(); ++itp)
        {
            number_t cr = itp->first;
            if ((at==_row && rc > cr) || (at==_col && rc < cr)) *(itv+(cr-1))= -conj(mat[itp->second]);
            else *(itv+(itp->first-1))=mat[itp->second];
        }
        break;
      case _selfAdjoint:
        for (;itp!=indadd.end(); ++itp)
        {
          number_t cr = itp->first;
          if ((at==_row && rc > cr) || (at==_col && rc < cr)) *(itv+(cr-1))= conj(mat[itp->second]);
          else *(itv+(itp->first-1))=mat[itp->second];
        }
        break;
      default :
        for (;itp!=indadd.end(); ++itp) *(itv+(itp->first-1))=mat[itp->second];
    }
    return v;
  }
  where("MatrixEntry::getRowCol(...)");
  error("scalar_only");
  return VectorEntry(_real,_scalar,1);  //fake return
}

FactorizationType MatrixEntry::factorization() const
{
   if(rEntries_p != 0)  return rEntries_p->factorization_;
   if(cEntries_p != 0)  return cEntries_p->factorization_;
   if(rmEntries_p != 0) return rmEntries_p->factorization_;
   if(cmEntries_p != 0) return cmEntries_p->factorization_;
   return _noFactorization;
}

// delete rows r1,..., r2 in matrix (may be expansive)
void MatrixEntry::deleteRows(number_t r1, number_t r2)
{
   if(rEntries_p != 0)  { rEntries_p->deleteRows(r1,r2);  return;}
   if(cEntries_p != 0)  { cEntries_p->deleteRows(r1,r2);  return;}
   if(rmEntries_p != 0) { rmEntries_p->deleteRows(r1,r2); return;}
   if(cmEntries_p != 0) { cmEntries_p->deleteRows(r1,r2); return;}
}

// delete cols c1,..., c2 in matrix (may be expansive)
void MatrixEntry::deleteCols(number_t c1, number_t c2)
{
   if(rEntries_p != 0)  { rEntries_p->deleteCols(c1,c2); return;}
   if(cEntries_p != 0)  { cEntries_p->deleteCols(c1,c2); return;}
   if(rmEntries_p != 0) { rmEntries_p->deleteCols(c1,c2); return;}
   if(cmEntries_p != 0) { cmEntries_p->deleteCols(c1,c2); return;}
}

//convert matrix storage to skyline storage
void MatrixEntry::toSkyline()
{
  if(rEntries_p != 0)  {rEntries_p->toSkyline(); return;}
  if(cEntries_p != 0)  {cEntries_p->toSkyline(); return;}
  if(rmEntries_p != 0) {rmEntries_p->toSkyline(); return;}
  if(cmEntries_p != 0) {cmEntries_p->toSkyline(); return;}
  where("MatrixEntry::toSkyline");
  error("null_pointer", "xxEntries_p");
}

// convert current matrix storage to an other matrix storage
void MatrixEntry::toStorage(MatrixStorage* ms)
{
  if(rEntries_p != 0)  {rEntries_p->toStorage(ms); return;}
  if(cEntries_p != 0)  {cEntries_p->toStorage(ms); return;}
  if(rmEntries_p != 0) {rmEntries_p->toStorage(ms); return;}
  if(cmEntries_p != 0) {cmEntries_p->toStorage(ms); return;}
  where("MatrixEntry::toStorage");
  error("null_pointer", "xxEntries_p");
}

// create new scalar matrix entry from non scalar matrix entry
MatrixEntry* MatrixEntry::toScalar()
{
   if(nbOfComponents == std::pair<dimen_t,dimen_t>(1,1)) return this;  //already scalar, nothing done

   //create new scalar entries
   MatrixEntry* mat=new MatrixEntry();
   mat->valueType_=valueType_;
   mat->strucType_=_scalar;
   mat->nbOfComponents= std::pair<dimen_t,dimen_t>(1,1);
   if(rmEntries_p != 0) {mat->rEntries_p=rmEntries_p->toScalar(real_t(0));}
   if(cmEntries_p != 0) {mat->cEntries_p=cmEntries_p->toScalar(complex_t(0));}
   return mat;
}

// convert current matrix to unsymmetric representation if it has a symmetric one
void MatrixEntry::toUnsymmetric(AccessType at)
{
  if(rEntries_p != 0)  {rEntries_p->toUnsymmetric(at); return;}
  if(cEntries_p != 0)  {cEntries_p->toUnsymmetric(at); return;}
  if(rmEntries_p != 0) {rmEntries_p->toUnsymmetric(at); return;}
  if(cmEntries_p != 0) {cmEntries_p->toUnsymmetric(at); return;}
  where("MatrixEntry::toUnsymmetric");
  error("null_pointer", "xxEntries_p");
}

// change real entries to complex entries : rEntries-> cEntries or rmEntries -> cmEntries
// nothing done when MatrixEntry is a complex one
MatrixEntry& MatrixEntry::toComplex()
{
   if(rEntries_p!=0)
   {
      cEntries_p= new LargeMatrix<complex_t>(rEntries_p->storagep(), complex_t(0), rEntries_p->sym);
      std::vector<real_t>::iterator itr=rEntries_p->values().begin();
      std::vector<complex_t>::iterator itc=cEntries_p->values().begin();
      for(; itr!=rEntries_p->values().end(); ++itr, ++itc) *itc= complex_t(*itr);
      delete rEntries_p;
      rEntries_p=0;
   }
   if(rmEntries_p!=0)
   {
      cmEntries_p= new LargeMatrix<Matrix<complex_t> >(rmEntries_p->storagep(), Matrix<complex_t>(0), rmEntries_p->sym);
      std::vector<Matrix<real_t> >::iterator itr=rmEntries_p->values().begin();
      std::vector<Matrix<complex_t> >::iterator itc=cmEntries_p->values().begin();
      for(; itr!=rmEntries_p->values().end(); ++itr, ++itc) *itc= cmplx(*itr);
      delete rmEntries_p;
      rmEntries_p=0;
   }
   valueType_=_complex;
   return *this;
}

// change complex entries to real entries : cEntries-> rEntries or cmEntries -> rmEntries
//   if realPart is true take real parts else take imag parts
// nothing done when MatrixEntry is a real one
MatrixEntry& MatrixEntry::toReal(bool realPart)
{
   if(cEntries_p!=0)
   {
      SymType sy = cEntries_p->sym;
      if(sy==_selfAdjoint) sy=_symmetric;
      if(sy==_skewAdjoint) sy=_skewSymmetric;
      rEntries_p= new LargeMatrix<real_t>(cEntries_p->storagep(), 0., sy);
      std::vector<complex_t>::iterator itc=cEntries_p->values().begin();
      std::vector<real_t>::iterator itr=rEntries_p->values().begin();
      if(realPart)
         for(; itc!=cEntries_p->values().end(); ++itr, ++itc) *itr= itc->real();
      else
         for(; itc!=cEntries_p->values().end(); ++itr, ++itc) *itr= itc->imag();
      delete cEntries_p;
      cEntries_p=0;
   }
   if(cmEntries_p!=0)
   {
      SymType sy = cmEntries_p->sym;
      if(sy==_selfAdjoint) sy=_symmetric;
      if(sy==_skewAdjoint) sy=_skewSymmetric;
      rmEntries_p= new LargeMatrix<Matrix<real_t> >(cmEntries_p->storagep(), Matrix<real_t>(0), sy);
      std::vector<Matrix<real_t> >::iterator itr=rmEntries_p->values().begin();
      std::vector<Matrix<complex_t> >::iterator itc=cmEntries_p->values().begin();
      if(realPart)
         for(; itc!=cmEntries_p->values().end(); ++itr, ++itc) *itr= itc->real();
      else
         for(; itc!=cmEntries_p->values().end(); ++itr, ++itc) *itr= itc->imag();
      delete cmEntries_p;
      cmEntries_p=0;
   }
   valueType_=_real;
   return *this;
}

// change entries to conjugate entries if entries are not real
// nothing done when MatrixEntry is a real one
MatrixEntry& MatrixEntry::toConj()
{
   if(cEntries_p!=0)  cEntries_p->toConj();
   if(cmEntries_p!=0) cmEntries_p->toConj();
   return *this;
}

// test if all non diagonal coefficients are zero
bool MatrixEntry::isDiagonal() const
{
   if(rEntries_p != 0)  return rEntries_p->isDiagonal();
   if(cEntries_p != 0)  return cEntries_p->isDiagonal();
   if(rmEntries_p != 0) return rmEntries_p->isDiagonal();
   if(cmEntries_p != 0) return cmEntries_p->isDiagonal();
   return true;
}

// test if all non diagonal coefficients are zero and diagonal coefficients are almost 1
bool MatrixEntry::isId(const double & tol) const
{
   if(rEntries_p != 0)  return rEntries_p->isId(tol);
   if(cEntries_p != 0)  return cEntries_p->isId(tol);
   if(rmEntries_p != 0) return rmEntries_p->isId(tol);
   if(cmEntries_p != 0) return cmEntries_p->isId(tol);
   return true;
}

// round to zero all coefficients close to 0 (|.| < aszero)
void MatrixEntry::roundToZero(real_t asZero)
{
   if(rEntries_p != 0) {  rEntries_p->roundToZero(asZero);  return; }
   if(cEntries_p != 0) {  cEntries_p->roundToZero(asZero);  return; }
   if(rmEntries_p != 0) { rmEntries_p->roundToZero(asZero); return; }
   if(cmEntries_p != 0) { cmEntries_p->roundToZero(asZero); return; }
}

//! set cols or rows to 0 (index start from 1)
void MatrixEntry::setColToZero(number_t c1, number_t c2)
{
   if(rEntries_p != 0)  rEntries_p->setColToZero(c1,c2);
   if(cEntries_p != 0)  cEntries_p->setColToZero(c1,c2);
   if(rmEntries_p != 0) rmEntries_p->setColToZero(c1,c2);
   if(cmEntries_p != 0) cmEntries_p->setColToZero(c1,c2);
}
void MatrixEntry::setRowToZero(number_t r1, number_t r2)
{
   if(rEntries_p != 0)  rEntries_p->setRowToZero(r1,r2);
   if(cEntries_p != 0)  cEntries_p->setRowToZero(r1,r2);
   if(rmEntries_p != 0) rmEntries_p->setRowToZero(r1,r2);
   if(cmEntries_p != 0) cmEntries_p->setRowToZero(r1,r2);
}

/*! change matrix entry representation : from current representation to matrix of nbr x nbc matrices
    rowPos : list of row indexes to be moved, nbr = rowPos.size
    colPos : list of col indexes to be moved, nbc = colPos.size
    for instance                                                                  [0 0]
        rowPos=[ 0 2 0], colPos=[1 2] means to move from representation [x y] to  [x y]
                                                                                  [0 0]
                                                                                  [0 0 x]
        rowPos=[ 1 0 0], colPos=[0 0 3] means to move from representation [x] to  [0 0 0]
                                                                                  [0 0 0]
    Note : rowPos and colPos has to be consistent with current representation (not checked here)

*/
void MatrixEntry::toMatrix(const std::vector<dimen_t>& rowPos, const std::vector<dimen_t>& colPos)
{
   dimen_t nbr=rowPos.size(), nbc=colPos.size();
   if(nbr==1 && nbc==1) return;     //nothing is done (no 1 x 1 matrices)
   dimPair dims(nbr, nbc);
   if(nbOfComponents==dims) return; //same representation

   if(rEntries_p!=0 || cEntries_p!=0)   //current is a matrix of real scalar
   {
      dimen_t r=nbr, c=nbc, k=0;
      while(k<nbr)  {if(rowPos[k]>0) {r=k; k=nbr;} else k++;}
      k=0;
      while(k<nbc)  {if(colPos[k]>0) {c=k; k=nbc;} else k++;}
      if(r==nbr || c==nbc) return;  //nothing is done, rowPos or colPos is null
      r++, c++;
      if(rEntries_p!=0)
      {
         Matrix<real_t> mat(nbr, nbc, 0.);
         rmEntries_p = new LargeMatrix<Matrix<real_t> >(rEntries_p->storagep(),mat,rEntries_p->sym);
         std::vector<real_t>::iterator it=rEntries_p->values().begin();
         std::vector<Matrix<real_t> >::iterator itm=rmEntries_p->values().begin();
         for(; it!=rEntries_p->values().end(); it++, itm++)
         {
            mat(r,c)=*it;
            *itm = mat;
         }
         delete rEntries_p;
         rEntries_p=0;
         strucType_=_matrix;
         nbOfComponents=dims;
         return;
      }
      else
      {
         Matrix<complex_t> mat(nbr, nbc, 0.);
         cmEntries_p = new LargeMatrix<Matrix<complex_t> >(cEntries_p->storagep(),mat,cEntries_p->sym);
         std::vector<complex_t>::iterator it=cEntries_p->values().begin();
         std::vector<Matrix<complex_t> >::iterator itm=cmEntries_p->values().begin();
         for(; it!=cEntries_p->values().end(); it++, itm++)
         {
            mat(r,c)=*it;
            *itm = mat;
         }
         delete cEntries_p;
         cEntries_p=0;
         strucType_=_matrix;
         nbOfComponents=dims;
         return;
      }
   }

   //current entry is matrix of matrices
   std::vector<dimen_t>::const_iterator itr, itc;
   if(rmEntries_p!=0)
   {
      Matrix<real_t> mat(nbr, nbc, 0.);
      std::vector<Matrix<real_t> >::iterator it=rmEntries_p->values().begin();
      for(; it!=rmEntries_p->values().end(); it++)
      {
         Matrix<real_t>::iterator itv=it->begin();
         for(dimen_t r=0; r<nbr; r++)
         {
            dimen_t i=rowPos[r];
            if(i!=0)
               for(dimen_t c=0; c<nbc; c++)
               {
                  dimen_t j=colPos[c];
                  if(j!=0) { mat(i,j)=*itv; itv++;}
               }
         }
         *it = mat;
      }
      nbOfComponents=dims;
      return;
   }
   if(cmEntries_p!=0)
   {
      Matrix<complex_t> mat(nbr, nbc, 0.);
      std::vector<Matrix<complex_t> >::iterator it=cmEntries_p->values().begin();
      for(; it!=cmEntries_p->values().end(); it++)
      {
         Matrix<complex_t>::iterator itv=it->begin();
         for(dimen_t r=0; r<nbr; r++)
         {
            dimen_t i=rowPos[r];
            if(i!=0)
               for(dimen_t c=0; c<nbc; c++)
               {
                  dimen_t j=colPos[c];
                  if(j!=0) { mat(i,j)=*itv; itv++;}
               }
         }
         *it = mat;
      }
      nbOfComponents=dims;
      return;
   }
}

/*!assign MatrixEntry m to current MatrixEntry
  num_r : row dof numbers of MatrixEntry m related to current MatrixEntry
  num_c : col dof numbers of MatrixEntry m related to current MatrixEntry
  if num_r (resp. num_c) is void, it means that the row (resp. col) dof numbers of MatrixEntry m is the same as current MatrixEntry
*/
void MatrixEntry::assign(const MatrixEntry& m, const std::vector<number_t>& num_r, const std::vector<number_t>& num_c)
{
  if(m.valueType_ == _real && valueType_ == _real) //real case
  {
    if(rEntries_p != 0)
    {
      if(m.rEntries_p != 0) {rEntries_p->assign(*m.rEntries_p, num_r, num_c); return;}
      where("MatrixEntry::assign");
      error("entry_inconsistent_structures");
    }
    if(rmEntries_p != 0)
    {
      if(m.rmEntries_p != 0) {rmEntries_p->assign(*m.rmEntries_p, num_r, num_c); return;}
      where("MatrixEntry::assign");
      error("entry_inconsistent_structures");
    }
    where("MatrixEntry::assign");
    error("null_pointer", "real Entries_p");
  }

  //complex case
  if(cEntries_p != 0)
  {
    if(m.rEntries_p != 0) {cEntries_p->assign(*m.rEntries_p, num_r, num_c); return;}
    if(m.cEntries_p != 0) {cEntries_p->assign(*m.cEntries_p, num_r, num_c); return;}
    where("MatrixEntry::assign");
    error("entry_inconsistent_structures");
  }
  if(cmEntries_p != 0)
  {
    if(m.rmEntries_p != 0) {cmEntries_p->assign(*m.rmEntries_p, num_r, num_c); return;}
    if(m.cmEntries_p != 0) {cmEntries_p->assign(*m.cmEntries_p, num_r, num_c); return;}
    where("MatrixEntry::assign");
    error("entry_inconsistent_structures");
  }
  where("MatrixEntry::assign");
  error("null_pointer", "complex Entries_p");
}


/*--------------------------------------------------------------------------------------------------------------------
   in/out operations
--------------------------------------------------------------------------------------------------------------------*/

// print entry
void MatrixEntry::viewStorage(std::ostream& out) const
{
   if(rEntries_p != 0)  rEntries_p->viewStorage(out);
   if(cEntries_p != 0)  cEntries_p->viewStorage(out);
   if(rmEntries_p != 0) rmEntries_p->viewStorage(out);
   if(cmEntries_p != 0) cmEntries_p->viewStorage(out);
}

// save matrix entries to file (dense or Matlab format)
void MatrixEntry::saveToFile(const string_t& fn, StorageType st, bool encode) const
{
   if(rEntries_p != 0)  rEntries_p->saveToFile(fn, st, encode);
   if(cEntries_p != 0)  cEntries_p->saveToFile(fn, st, encode);
   if(rmEntries_p != 0) rmEntries_p->saveToFile(fn, st, encode);
   if(cmEntries_p != 0) cmEntries_p->saveToFile(fn, st, encode);
}

/*--------------------------------------------------------------------------------------------------------------------
   algebraic operations
--------------------------------------------------------------------------------------------------------------------*/
MatrixEntry& MatrixEntry::operator*=(const real_t& r)
{
   if(rEntries_p != 0)  *rEntries_p *= r;
   if(cEntries_p != 0)  *cEntries_p *= r;
   if(rmEntries_p != 0) *rmEntries_p *= r;
   if(cmEntries_p != 0) *cmEntries_p *= r;
   return *this;
}

MatrixEntry& MatrixEntry::operator*=(const complex_t& c)
{
   if(rEntries_p != 0)
   {
      cEntries_p = new LargeMatrix<complex_t>(*rEntries_p);
      delete rEntries_p;
      rEntries_p=0;
   }
   if(cEntries_p != 0)  *cEntries_p *= c;
   if(rmEntries_p != 0)
   {
      cmEntries_p = new LargeMatrix<Matrix<complex_t> >(*rmEntries_p);
      delete rmEntries_p;
      rmEntries_p=0;
   }
   if(cmEntries_p != 0) *cmEntries_p *= c;
   valueType_=_complex;
   return *this;
}

MatrixEntry& MatrixEntry::operator/=(const real_t& r)
{
   return (*this) *= (1. / r);
}

MatrixEntry& MatrixEntry::operator/=(const complex_t& c)
{
   return (*this) *= (1. / c);
}

//+= and -= operation for same storage ans same structure !!!
MatrixEntry& MatrixEntry::operator+=(const MatrixEntry& m)
{
  if(strucType_ != m.strucType_)
  {
    where("MatrixEntry::operator+=");
    error("entry_mismatch_structures", words("structure",strucType_), words("structure", m.strucType_));
  }
  if(storagep() != m.storagep())
  {
    where("MatrixEntry::operator+=");
    error("matrixentry_mismatch_storages");
  }

  if(rEntries_p != 0) //current is a scalar matrix
  {
    if(m.rEntries_p != 0) {*rEntries_p += *m.rEntries_p; return *this;}
    if(m.cEntries_p != 0) //add a complex matrix to a real one
    {
      //copy rEntries_p to cEntries_p
      cEntries_p = new LargeMatrix<complex_t>(*rEntries_p);
      delete rEntries_p;
      rEntries_p=0;
    }
  }
  if(cEntries_p != 0 && m.cEntries_p != 0) {*cEntries_p += *m.cEntries_p; return *this;}
  if(rmEntries_p != 0) //current is a scalar matrix
  {
    if(m.rmEntries_p != 0) {*rmEntries_p += *m.rmEntries_p; return *this;}
    if(m.cmEntries_p != 0) //add a complex matrix to a real one
    {
      //copy rEntries_p to cEntries_p
      cmEntries_p = new LargeMatrix<Matrix<complex_t> >(*rmEntries_p);
      delete rmEntries_p;
      rmEntries_p=0;
    }
  }
  if(cmEntries_p != 0 && m.cmEntries_p != 0) {*cmEntries_p += *m.cmEntries_p;}
  return *this;
}

MatrixEntry& MatrixEntry::operator-=(const MatrixEntry& m)
{
  if(strucType_ != m.strucType_)
  {
    where("MatrixEntry::operator-=");
    error("entry_mismatch_structures", words("structure",strucType_), words("structure", m.strucType_));
  }
  if(storagep() != m.storagep())
  {
    where("MatrixEntry::operator-=");
    error("matrixentry_mismatch_storages");
  }

  if(rEntries_p != 0) //current is a scalar matrix
  {
    if(m.rEntries_p != 0) {*rEntries_p -= *m.rEntries_p; return *this;}
    if(m.cEntries_p != 0) //add a complex matrix to a real one
    {
      //copy rEntries_p to cEntries_p
      cEntries_p = new LargeMatrix<complex_t>(*rEntries_p);
      delete rEntries_p;
      rEntries_p=0;
    }
  }
  if(cEntries_p != 0 && m.cEntries_p != 0) {*cEntries_p -= *m.cEntries_p; return *this;}
  if(rmEntries_p != 0) //current is a scalar matrix
  {
    if(m.rmEntries_p != 0) {*rmEntries_p -= *m.rmEntries_p; return *this;}
    if(m.cmEntries_p != 0) //add a complex matrix to a real one
    {
      //copy rEntries_p to cEntries_p
      cmEntries_p = new LargeMatrix<Matrix<complex_t> >(*rmEntries_p);
      delete rmEntries_p;
      rmEntries_p=0;
    }
  }
  if(cmEntries_p != 0 && m.cmEntries_p != 0) {*cmEntries_p -= *m.cmEntries_p;}
  return *this;
}

/*!
  add MatrixEntry m to current MatrixEntry with scalar multiplication
  num_r : row dof numbers of MatrixEntry m related to current MatrixEntry
  num_c : col dof numbers of MatrixEntry m related to current MatrixEntry
  if num_r (resp. num_c) is void, it means that the row (resp. col) dof numbers of MatrixEntry m is the same as current MatrixEntry
*/

void MatrixEntry::add(const MatrixEntry& m, const std::vector<number_t>& num_r, const std::vector<number_t>& num_c, complex_t a)
{
  if(a.imag() == 0 && m.valueType_ == _real && valueType_ == _real) //real case
  {
    if(rEntries_p != 0)
    {
      if(m.rEntries_p != 0) {rEntries_p->add(*m.rEntries_p, num_r, num_c, a.real()); return;}
      where("MatrixEntry::add");
      error("entry_inconsistent_structures");
    }
    if(rmEntries_p != 0)
    {
      if(m.rmEntries_p != 0) {rmEntries_p->add(*m.rmEntries_p, num_r, num_c, a.real()); return;}
      where("MatrixEntry::add");
      error("entry_inconsistent_structures");
    }
    where("MatrixEntry::add");
    error("null_pointer", "real Entries_p");
  }

  //complex case
  if(cEntries_p != 0)
  {
    if(m.rEntries_p != 0) {cEntries_p->add(*m.rEntries_p, num_r, num_c, a); return;}
    if(m.cEntries_p != 0) {cEntries_p->add(*m.cEntries_p, num_r, num_c, a); return;}
    where("MatrixEntry::add");
    error("entry_inconsistent_structures");
  }
  if(cmEntries_p != 0)
  {
    if(m.rmEntries_p != 0) {cmEntries_p->add(*m.rmEntries_p, num_r, num_c, a); return;}
    if(m.cmEntries_p != 0) {cmEntries_p->add(*m.cmEntries_p, num_r, num_c, a); return;}
    where("MatrixEntry::add");
    error("entry_inconsistent_structures");
  }
  where("MatrixEntry::add");
  error("null_pointer", "complex Entries_p");
}

//add entries (itb to ite) to current MatrixEntry (has to be consistent, no check)
template<typename iterator>
void MatrixEntry::add(iterator& itb, iterator& ite)
{
   if(rEntries_p != 0)  {rEntries_p->add(itb,ite); return;}
   if(rmEntries_p != 0) {rmEntries_p->add(itb,ite); return;}
   if(cEntries_p != 0)  {cEntries_p->add(itb,ite); return;}
   if(cmEntries_p != 0) {cmEntries_p->add(itb,ite); return;}
}

/*! combine cols of matrix : c2 = c2 + a *c1 (indices start from 1)
    if updatestorage is true, the storage is dynamically updated (very expansive operation !!!)
*/
void MatrixEntry::addColToCol(number_t c1, number_t c2, complex_t a, bool updateStorage)
{
   if(a==complex_t(0)) return;  //nothing to add
   if(c1<1 || c1>nbOfCols() || c2<1 || c2>nbOfCols()) return; //index out of bounds

   if(rEntries_p != 0)  {rEntries_p->addColToCol(c1,c2,a.real(), updateStorage); return;}
   if(rmEntries_p != 0) {rmEntries_p->addColToCol(c1,c2,a.real(), updateStorage); return;}
   if(cEntries_p != 0)  {cEntries_p->addColToCol(c1,c2,a, updateStorage); return;}
   if(cmEntries_p != 0) {cmEntries_p->addColToCol(c1,c2,a, updateStorage); return;}
}

/*! combine rows of matrix : r2 = rc2 + a *r1 (indices start from 1)
    if updatestorage is true, the storage is dynamically updated (very expansive operation !!!)
*/
void MatrixEntry::addRowToRow(number_t r1, number_t r2, complex_t a, bool updateStorage)
{
   if(a==complex_t(0)) return;  //nothing to add
   if(r1<1 || r1>nbOfRows() || r2<1 || r2>nbOfRows()) return; //index out of bounds

   if(rEntries_p != 0)  {rEntries_p->addRowToRow(r1,r2,a.real(), updateStorage); return;}
   if(rmEntries_p != 0) {rmEntries_p->addRowToRow(r1,r2,a.real(), updateStorage); return;}
   if(cEntries_p != 0)  {cEntries_p->addRowToRow(r1,r2,a, updateStorage); return;}
   if(cmEntries_p != 0) {cmEntries_p->addRowToRow(r1,r2,a, updateStorage); return;}
}

/*! copy in current matrix values of mat located at mat_adrs at position cur_adrs
    it is assumed that mat and current matrix are consistent
        in structure type scalar/scalar or matrix/matrix
        in value type real/real  or complex/complex
    if not do conversion before
    it is also assumed that mat_adrs and cur_adrs have the same size (addresses map!)
*/
void MatrixEntry::copyVal(const MatrixEntry& mat, const std::vector<number_t>& mat_adrs,const std::vector<number_t>& cur_adrs)
{
  if (rEntries_p != 0 && mat.rEntries_p!=0) {rEntries_p->copyVal(*mat.rEntries_p, mat_adrs, cur_adrs); return;}
  if (cEntries_p != 0 && mat.cEntries_p!=0) {cEntries_p->copyVal(*mat.cEntries_p, mat_adrs, cur_adrs); return;}
  if (rmEntries_p != 0 && mat.rmEntries_p!=0) {rmEntries_p->copyVal(*mat.rmEntries_p, mat_adrs, cur_adrs); return;}
  if (cmEntries_p != 0 && mat.cmEntries_p!=0) {cmEntries_p->copyVal(*mat.cmEntries_p, mat_adrs, cur_adrs); return;}
  where("MatrixEntry::copyVal");
  error("entry_inconsistent_structures");
}


/*! partial norm of column c, only components in range [r1,r2], say
    sqrt( sum_{r1 <= r <= r2} Mrc^2 )
*/
real_t MatrixEntry::partialNormOfCol(number_t c, number_t r1, number_t r2) const
{
   if(rEntries_p != 0)  return rEntries_p->partialNormOfCol(c,r1,r2);
   if(rmEntries_p != 0) return rmEntries_p->partialNormOfCol(c,r1,r2);
   if(cEntries_p != 0)  return cEntries_p->partialNormOfCol(c,r1,r2);
   if(cmEntries_p != 0) return cmEntries_p->partialNormOfCol(c,r1,r2);
   return 0.;
}

// Frobenius norm
real_t MatrixEntry::norm2() const
{
   if(rEntries_p != 0)  return rEntries_p->norm2();
   if(rmEntries_p != 0) return rmEntries_p->norm2();
   if(cEntries_p != 0)  return cEntries_p->norm2();
   if(cmEntries_p != 0) return cmEntries_p->norm2();
   return 0.;
}

// infinite norm
real_t MatrixEntry::norminfty() const
{
   if(rEntries_p != 0)  return rEntries_p->norminfty();
   if(rmEntries_p != 0) return rmEntries_p->norminfty();
   if(cEntries_p != 0)  return cEntries_p->norminfty();
   if(cmEntries_p != 0) return cmEntries_p->norminfty();
   return 0.;
}

// product of MatrixEntry
MatrixEntry operator*(const MatrixEntry& mA,const MatrixEntry& mB)
{
  if(mA.nbOfCols()!= mB.nbOfRows() || mA.nbOfComponents.second != mB.nbOfComponents.first)
  {
    where("MatrixEntry * MatrixEntry");
    error("entry_inconsistent_structures");
  }
  ValueType vt=_real;
  if(mA.valueType_==_complex || mB.valueType_==_complex) vt=_complex;
  StrucType st=_scalar;
  dimen_t nbrs=mA.nbOfComponents.first, nbcs=mB.nbOfComponents.second;
  if(nbrs!=1 || nbcs!=1) st=_matrix;
  MatrixEntry mR(vt,st,0,dimPair(nbrs,nbcs),_noSymmetry);

  switch(st)
  {
    case _scalar :
      switch(vt)
    {
      case _real :    mR.rEntries_p = new LargeMatrix<real_t>(0,_noSymmetry);
        multMatrixMatrix<real_t, real_t, real_t>(*mA.rEntries_p,*mB.rEntries_p, *mR.rEntries_p);
        return mR;
      case _complex :
      {
        if(mA.rEntries_p!=0)
        {
          mR.cEntries_p = new LargeMatrix<complex_t>(0,_noSymmetry);
          multMatrixMatrix<real_t, complex_t, complex_t>(*mA.rEntries_p, *mB.cEntries_p, *mR.cEntries_p);
          return mR;
        }
        if(mB.rEntries_p!=0)
        {
          mR.cEntries_p = new LargeMatrix<complex_t>(0,_noSymmetry);
          multMatrixMatrix<complex_t, real_t, complex_t>(*mA.cEntries_p, *mB.rEntries_p, *mR.cEntries_p);
          return mR;
        }
        else
        {
          mR.cEntries_p = new LargeMatrix<complex_t>(0,_noSymmetry);
          multMatrixMatrix<complex_t,complex_t, complex_t>(*mA.cEntries_p, *mB.cEntries_p, *mR.cEntries_p);
          return mR;
        }
      }
      default : error("matrixentry_abnormal_type", words("value", vt));
    }
//      case _matrix :
//        switch(vt)
//          {
//            case _real :
//              {
//                mR.rmEntries_p = new LargeMatrix<Matrix<real_t> >(0, dimPair(nbrs,nbcs),_noSymmetry);
//                multMatrixMatrix<Matrix<real_t>,Matrix<real_t>,Matrix<real_t> >(*mA.rmEntries_p, *mB.rmEntries_p, *mR.rmEntries_p);
//                return mR;
//              }
//            case _complex :
//              {
//                if(mA.rmEntries_p!=0)
//                  {
//                    mR.cmEntries_p = new LargeMatrix<Matrix<complex_t> >(0, dimPair(nbrs,nbcs),_noSymmetry);
//                    multMatrixMatrix<Matrix<real_t>,Matrix<complex_t>,Matrix<complex_t> >(*mA.rmEntries_p, *mB.cmEntries_p, *mR.cmEntries_p);
//                    return mR;
//                  }
//                if(mB.rmEntries_p!=0)
//                  {
//                    mR.cmEntries_p = new LargeMatrix<Matrix<complex_t> >(0, dimPair(nbrs,nbcs),_noSymmetry);
//                    multMatrixMatrix<Matrix<complex_t>,Matrix<real_t>,Matrix<complex_t> >(*mA.cmEntries_p, *mB.rmEntries_p, *mR.cmEntries_p);
//                    return mR;
//                  }
//                else
//                  {
//                    mR.cmEntries_p = new LargeMatrix<Matrix<complex_t> >(0, dimPair(nbrs,nbcs),_noSymmetry);
//                    multMatrixMatrix<Matrix<complex_t>,Matrix<complex_t>,Matrix<complex_t> >(*mA.cmEntries_p, *mB.cmEntries_p, *mR.cmEntries_p);
//                    return mR;
//                  }
//              }
//            default : error("matrixentry_abnormal_type", words("value", vt));
//          }
      default : error("matrixentry_abnormal_type", words("structure", st));
   }

   return mR;
}

//! return the conjugate of a matrix entry
MatrixEntry conj(const MatrixEntry & me)
{
   if(me.rEntries_p != 0 || me.rmEntries_p != 0)  {return me;}
   MatrixEntry mecj=me;
   if(mecj.cEntries_p!=0) mecj.cEntries_p->toConj();
   if(mecj.cmEntries_p!=0) mecj.cmEntries_p->toConj();
   return mecj;
}

/*--------------------------------------------------------------------------------------------------------------------
   factorization operations
--------------------------------------------------------------------------------------------------------------------*/
//factorize matrix entry as LU,
void MatrixEntry::luFactorize(bool withPermutation)
{
  if(rEntries_p != 0)  {rEntries_p->luFactorize(withPermutation); return;}
  if(cEntries_p != 0)  {cEntries_p->luFactorize(withPermutation); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::luFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::luFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::luFactorize");
  error("null_pointer", "xxEntries_p");
}
//factorize matrix entry as LDLt,
void MatrixEntry::ldltFactorize()
{
  if(rEntries_p != 0)  {rEntries_p->ldltFactorize(); return;}
  if(cEntries_p != 0)  {cEntries_p->ldltFactorize(); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::ldltFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::ldltFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::ldltFactorize");
  error("null_pointer", "xxEntries_p");
}

//factorize matrix entry as LDLt,
void MatrixEntry::ldlstarFactorize()
{
  if(rEntries_p != 0)  {rEntries_p->ldlstarFactorize(); return;}
  if(cEntries_p != 0)  {cEntries_p->ldlstarFactorize(); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::ldlstarFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::ldlstarFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::ldlstarFactorize");
  error("null_pointer", "xxEntries_p");
}


//incomplete factorization operations
//incomplete factorize matrix entry as LU,
void MatrixEntry::iluFactorize()
{
  if(rEntries_p != 0)  {rEntries_p->iluFactorize(); return;}
  if(cEntries_p != 0)  {cEntries_p->iluFactorize(); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::iluFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::iluFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::iluFactorize");
  error("null_pointer", "xxEntries_p");
}

//incomplete factorize matrix entry as LDLt,
void MatrixEntry::ildltFactorize()
{
  if(rEntries_p != 0)  {rEntries_p->ildltFactorize(); return;}
  if(cEntries_p != 0)  {cEntries_p->ildltFactorize(); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::ildltFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::ildltFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::ildltFactorize");
  error("null_pointer", "xxEntries_p");
}

//incomplete factorize matrix entry as LDL*,
void MatrixEntry::ildlstarFactorize()
{
  if(rEntries_p != 0)  {rEntries_p->ildlstarFactorize(); return;}
  if(cEntries_p != 0)  {cEntries_p->ildlstarFactorize(); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::ildlstarFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::ildlstarFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::ildlstarFactorize");
  error("null_pointer", "xxEntries_p");
}

//incomplete factorize matrix entry as LLt,
void MatrixEntry::illtFactorize()
{
  if(rEntries_p != 0)  {rEntries_p->illtFactorize(); return;}
  if(cEntries_p != 0)  {cEntries_p->illtFactorize(); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::illtFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::illtFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::illtFactorize");
  error("null_pointer", "xxEntries_p");
}

//incomplete factorize matrix entry as LL*,
void MatrixEntry::illstarFactorize()
{
  if(rEntries_p != 0)  {rEntries_p->illstarFactorize(); return;}
  if(cEntries_p != 0)  {cEntries_p->illstarFactorize(); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::illstarFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::illstarFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::illstarFactorize");
  error("null_pointer", "xxEntries_p");
}
//factorize matrix using umfpack
void MatrixEntry::umfpackFactorize()
{
  if(rEntries_p != 0)  {rEntries_p->umfpackFactorize(); return;}
  if(cEntries_p != 0)  {cEntries_p->umfpackFactorize(); return;}
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::umfpackFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::umfpackFactorize");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::umfpackFactorize");
  error("null_pointer", "xxEntries_p");
}

//solve linear system using LDLt factorization
void MatrixEntry::ldltSolve(const VectorEntry& B, VectorEntry& X)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->ldltSolve(*B.rEntries_p, *X.rEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->ldltSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::ldltSolve");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->ldltSolve(*B.rEntries_p, *X.cEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->ldltSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::ldltSolve");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::ldltSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::ldltSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::ldltSolve");
  error("null_pointer", "xxEntries_p");
}

//solve linear system using LDL* factorization
void MatrixEntry::ldlstarSolve(const VectorEntry& B, VectorEntry& X)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->ldlstarSolve(*B.rEntries_p, *X.rEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->ldlstarSolve(*B.rEntries_p, *X.rEntries_p); return;}
    where("MatrixEntry::ldlstarSolve");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->ldlstarSolve(*B.rEntries_p, *X.cEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->ldlstarSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::ldlstarSolve");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::ldlstarSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::ldlstarSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::ldlstarSolve");
  error("null_pointer", "xxEntries_p");
}

//solve linear system using LLt factorization
void MatrixEntry::lltSolve(const VectorEntry& B, VectorEntry& X)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->lltSolve(*B.rEntries_p, *X.rEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->lltSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::ldltSolve");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->lltSolve(*B.rEntries_p, *X.cEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->lltSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::lltSolve");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::lltSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::lltSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::lltSolve");
  error("null_pointer", "xxEntries_p");
}

//solve linear system using LU factorization
void MatrixEntry::luSolve(const VectorEntry& B, VectorEntry& X)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->luSolve(*B.rEntries_p, *X.rEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->luSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::luSolve");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->luSolve(*B.rEntries_p, *X.cEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->luSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::luSolve");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::luSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::luSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::luSolve");
  error("null_pointer", "xxEntries_p");
}

//solve linear system using umfpack LU factorization
void MatrixEntry::umfluSolve(const VectorEntry& B, VectorEntry& X)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->umfluSolve(*B.rEntries_p, *X.rEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->umfluSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::umfluSolve");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->umfluSolve(*B.rEntries_p, *X.cEntries_p); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->umfluSolve(*B.cEntries_p, *X.cEntries_p); return;}
    where("MatrixEntry::umfluSolve");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::umfluSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::umfluSolve");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::umfluSolve");
  error("null_pointer", "xxEntries_p");
}

//solve linear system using Sor Diagonal Solver
void MatrixEntry::sorDiagonalSolve(const VectorEntry& B,  VectorEntry& X, const real_t w)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->sorDiagonalSolve(*B.rEntries_p, *X.rEntries_p, w); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->sorDiagonalSolve(*B.cEntries_p, *X.cEntries_p, w); return;}
    where("MatrixEntry::sorDiagonalSolver");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->sorDiagonalSolve(*B.rEntries_p, *X.cEntries_p, w); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->sorDiagonalSolve(*B.cEntries_p, *X.cEntries_p, w); return;}
    where("MatrixEntry::sorDiagonalSolver");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::sorDiagonalSolver");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::sorDiagonalSolver");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::sorDiagonalSolver");
  error("null_pointer", "xxEntries_p");
}
//solve linear system using Sor Upper Solver
void MatrixEntry::sorDiagonalMatrixVector(const VectorEntry& B,  VectorEntry& X, const real_t w)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->sorDiagonalMatrixVector(*B.rEntries_p, *X.rEntries_p, w); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->sorDiagonalMatrixVector(*B.cEntries_p, *X.cEntries_p, w); return;}
    where("MatrixEntry::sorDiagonalMatrixVector");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->sorDiagonalMatrixVector(*B.rEntries_p, *X.cEntries_p, w); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->sorDiagonalMatrixVector(*B.cEntries_p, *X.cEntries_p, w); return;}
    where("MatrixEntry::sorDiagonalMatrixVector");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::sorDiagonalMatrixVector");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::sorDiagonalMatrixVector");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::sorDiagonalMatrixVector");
  error("null_pointer", "xxEntries_p");
}

//solve linear system using Sor Upper Solver
void MatrixEntry::sorUpperSolve(const VectorEntry& B,  VectorEntry& X, const real_t w)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->sorUpperSolve(*B.rEntries_p, *X.rEntries_p, w); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->sorUpperSolve(*B.cEntries_p, *X.cEntries_p, w); return;}
    where("MatrixEntry::sorUpperSolver");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->sorUpperSolve(*B.rEntries_p, *X.cEntries_p, w); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->sorUpperSolve(*B.cEntries_p, *X.cEntries_p, w); return;}
    where("MatrixEntry::sorUpperSolver");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::sorUpperSolver");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::sorUpperSolver");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::sorUpperSolver");
  error("null_pointer", "xxEntries_p");
}

//solve linear system using Sor Upper Solver
void MatrixEntry::sorLowerSolve(const VectorEntry& B,  VectorEntry& X, const real_t w)
{
  if(rEntries_p != 0)  //scalar real matrix
  {
    if(B.rEntries_p != 0 && X.rEntries_p != 0) {rEntries_p->sorLowerSolve(*B.rEntries_p, *X.rEntries_p, w); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {rEntries_p->sorLowerSolve(*B.cEntries_p, *X.cEntries_p, w); return;}
    where("MatrixEntry::sorLowerSolver");
    error("entry_inconsistent_structures");
  }
  if(cEntries_p != 0)  //scalar complex matrix
  {
    if(B.rEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->sorLowerSolve(*B.rEntries_p, *X.cEntries_p, w); return;}
    if(B.cEntries_p != 0 && X.cEntries_p != 0) {cEntries_p->sorLowerSolve(*B.cEntries_p, *X.cEntries_p, w); return;}
    where("MatrixEntry::sorLowerSolver");
    error("entry_inconsistent_structures");
  }
  if(rmEntries_p != 0)
  {
    where("MatrixEntry::sorLowerSolver");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(cmEntries_p != 0)
  {
    where("MatrixEntry::sorLowerSolver");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  where("MatrixEntry::sorLowerSolver");
  error("null_pointer", "xxEntries_p");
}

std::ostream& operator<<(std::ostream& out, const MatrixEntry& mat)
{
  mat.print(out);
  return out;
}

//---------------------------------------------------------------------------------------------
//! matrix * vector (consistent structure)
//---------------------------------------------------------------------------------------------
VectorEntry operator*(const MatrixEntry& mat, const VectorEntry& vec)
{
  VectorEntry res;
  multMatrixVector(mat,vec,res);
  return res;
}

void multMatrixVector(const MatrixEntry& mat, const VectorEntry& vec, VectorEntry& res)
{
  StrucType stm=mat.strucType_, stv=vec.strucType_;
  if ((stm==_scalar && stv!=_scalar) || (stm==_matrix && stv!=_vector))
  {
    where("MatrixEntry * VectorEntry");
    error("entry_inconsistent_structures");
  }
  number_t nv=1;
  if (stm==_matrix)
  {
    nv=mat.nbOfComponents.first;
    if(mat.nbOfComponents.second!=vec.nbOfComponents)
    {
      where("MatrixEntry * VectorEntry");
      error("entry_mismatch_dims", mat.nbOfComponents.second, vec.nbOfComponents);
    }
  }
  if (mat.nbOfCols()!=vec.size())
  {
    where("MatrixEntry * VectorEntry");
    error("entry_mismatch_dims", mat.nbOfCols(), vec.size());
  }
  ValueType vtm=mat.valueType_, vtv=vec.valueType_, vtr=_real;
  if (vtm==_complex || vtv==_complex) vtr=_complex;

  //reallocate res if not suited
  if (res.valueType_!=vtr || res.strucType_!=stm || res.size()!=mat.nbOfRows() ||res.nbOfComponents !=nv )
     res=VectorEntry(vtr, stm, mat.nbOfRows(),nv);

  if (stm==_scalar)
  {
    if (vtm==_real && vtv==_real)       *res.rEntries_p = *mat.rEntries_p * *vec.rEntries_p;
    if (vtm==_complex && vtv==_complex) *res.cEntries_p = *mat.cEntries_p * *vec.cEntries_p;
    if (vtm==_complex && vtv==_real)    *res.cEntries_p = *mat.cEntries_p * cmplx(*vec.rEntries_p);
    if (vtm==_real && vtv==_complex)
    {
      *res.cEntries_p = Vector<complex_t>(*mat.rEntries_p * real(*vec.cEntries_p));
      *res.cEntries_p += complex_t(0.,1.)* Vector<complex_t>(*mat.rEntries_p * imag(*vec.cEntries_p));
    }
  }
  else
  {
    if (vtm==_real && vtv==_real) *res.rvEntries_p = *mat.rmEntries_p * *vec.rvEntries_p;
    if (vtm==_complex && vtv==_complex) *res.cvEntries_p = *mat.cmEntries_p * *vec.cvEntries_p;
    if (vtm==_complex && vtv==_real) *res.cvEntries_p = *mat.cmEntries_p * cmplx(*vec.rvEntries_p);
    if (vtm==_real && vtv==_complex)
    {
      *res.cvEntries_p =  Vector<Vector<complex_t> >(*mat.rmEntries_p * real(*vec.cvEntries_p));
      *res.cvEntries_p += complex_t(0.,1.)* Vector<Vector<complex_t> >(*mat.rmEntries_p * imag(*vec.cvEntries_p));
    }
  }
}

//---------------------------------------------------------------------------------------------
//! vector * matrix (consistent structure)
//---------------------------------------------------------------------------------------------
VectorEntry operator*(const VectorEntry& vec, const MatrixEntry& mat)
{
  VectorEntry res;
  multVectorMatrix(vec,mat,res);
  return res;
}

void multVectorMatrix(const MatrixEntry& mat, const VectorEntry& vec, VectorEntry& res)
{
  multVectorMatrix(vec, mat, res);
}

void multVectorMatrix(const VectorEntry& vec, const MatrixEntry& mat, VectorEntry& res)
{
  StrucType stm=mat.strucType_, stv=vec.strucType_;
  if ((stm==_scalar && stv!=_scalar) || (stm==_matrix && stv!=_vector))
  {
    where("VectorEntry * MatrixEntry");
    error("entry_inconsistent_structures");
  }
  number_t nu=1;
  if (stm==_matrix)
  {
    nu=mat.nbOfComponents.second;
    if (mat.nbOfComponents.first!=vec.nbOfComponents)
    {
      where("VectorEntry * MatrixEntry");
      error("entry_mismatch_dims", vec.nbOfComponents, mat.nbOfComponents.first);
    }
  }
  if (mat.nbOfRows()!=vec.size())
  {
    where("VectorEntry * MatrixEntry");
    error("entry_mismatch_dims", vec.size(), mat.nbOfRows());
  }
  ValueType vtm=mat.valueType_, vtv=vec.valueType_, vtr=_real;
  if (vtm==_complex || vtv==_complex) vtr=_complex;

  //reallocate res if not suited
  if (res.valueType_!=vtr || res.strucType_!=stm || res.size()!=mat.nbOfCols() ||res.nbOfComponents !=nu )
     res=VectorEntry(vtr, stm, mat.nbOfCols(),nu);

  if (stm==_scalar)
  {
    if (vtm==_real && vtv==_real)       *res.rEntries_p = *vec.rEntries_p * *mat.rEntries_p;
    if (vtm==_complex && vtv==_complex) *res.cEntries_p = *vec.cEntries_p * *mat.cEntries_p;
    if (vtm==_complex && vtv==_real)    *res.cEntries_p = cmplx(*vec.rEntries_p) * *mat.cEntries_p;
    if (vtm==_real && vtv==_complex)
    {
      *res.cEntries_p = Vector<complex_t>(real(*vec.cEntries_p) * *mat.rEntries_p);
      *res.cEntries_p += complex_t(0.,1.)* Vector<complex_t>(imag(*vec.cEntries_p) * *mat.rEntries_p);
    }
  }
  else
  {
    if (vtm==_real && vtv==_real) *res.rvEntries_p = *vec.rvEntries_p * *mat.rmEntries_p;
    if (vtm==_complex && vtv==_complex) *res.cvEntries_p = *vec.cvEntries_p * *mat.cmEntries_p;
    if (vtm==_complex && vtv==_real) *res.cvEntries_p = cmplx(*vec.rvEntries_p) * *mat.cmEntries_p;
    if (vtm==_real && vtv==_complex)
    {
      *res.cvEntries_p =  Vector<Vector<complex_t> >(real(*vec.cvEntries_p) * *mat.rmEntries_p);
      *res.cvEntries_p += complex_t(0.,1.)* Vector<Vector<complex_t> >(imag(*vec.cvEntries_p) * *mat.rmEntries_p);
    }
  }
}


//---------------------------------------------------------------------------------------------
/*!QR factorization of rectangular matrix pxn with p<=n 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
   matQ     : unitary matrix of size p x p  stored as RowCsStorage
   computeQ : if true the matrix Q is computed else not
   rhs      : right hand side vector pointer (if 0 not reduction for rhs)
   withColPermutation : if true, column may be permuted (more stable factorization), may be forced by algorithm
   numcol_p : renumbering of columns if they have been permuted (pointer)
   stop     : iteration number when the algorithm stops: stop = nbr when algorithm exits normally */
//---------------------------------------------------------------------------------------------
void QR(const MatrixEntry& mat, MatrixEntry& matR, bool computeQ, MatrixEntry& matQ, VectorEntry* rhs,
        bool& withColPermutation, std::vector<number_t>*& numcol, number_t& stop)
{
  //initialization
  matR.clear();
  matQ.clear();
  ValueType vtm=mat.valueType_; StrucType stm=mat.strucType_; dimPair nbc=mat.nbOfComponents;
  matR.valueType_ = matQ.valueType_ = vtm;
  matR.strucType_ = matQ.strucType_ = stm;
  matR.nbOfComponents = matQ.nbOfComponents = nbc;
  if(rhs!=0)
  {
    if(rhs->strucType_ != _scalar)
    {
      where("QR(MatrixEntry, ...)");
      error("scalar_rhs_only");
    }

  }

  std::vector<real_t>* rsmb=0;
  std::vector<complex_t>* csmb=0;

  //call QR(LargeMatrix,...)
  if(stm == _scalar)
  {
    switch(vtm)
    {
      case _real :
      {
        if(rhs!= 0)
        {
            if(rhs->valueType_ == _real)  QR(*mat.rEntries_p, matR.rEntries_p, computeQ, matQ.rEntries_p, rhs->rEntries_p,
                                             withColPermutation, numcol, stop);   //QR factorisation using Housholder method
            else                         QR(*mat.rEntries_p, matR.rEntries_p, computeQ, matQ.rEntries_p, rhs->cEntries_p,
                                            withColPermutation, numcol, stop);   //QR factorisation using Housholder method
        }
        else QR(*mat.rEntries_p, matR.rEntries_p, computeQ, matQ.rEntries_p, rsmb,
                withColPermutation, numcol, stop);   //QR factorisation using Housholder method
      }
      break;
      default : //complex case
      {
        if(rhs!=0 && rhs->valueType_ == _real) rhs->toComplex(); //convert rhs to complex
        if(rhs!=0) csmb = rhs->cEntries_p;
        QR(*mat.cEntries_p, matR.cEntries_p, computeQ, matQ.cEntries_p, csmb,
           withColPermutation, numcol, stop);   //QR factorisation using Housholder method
      }
    }
  }
  else
  {
    where("QR(MatrixEntry, ...)");
    error("matrixentry_matrixofmatrices_not_handled");
  }
}

//-----------------------------------------------------------------------------------------------------------------------------
/*! solve an pxp upper triangular system
    mat      : pxp upper triangular matrix with its diagonal coeffs = 1,
    matR     : rectangular matrix p x m, understood as a list of right hand side vectors of size p (if 0 no computation for matR)
    rhs      : a unique right hand side vector pointer (if 0 no computation for rhs)

    At the end, rhs contains the vector x solution of mat * x = rhs and matR contains the vectors x1, x2, ...xm solutions of mat * xj = matRj
    Note: at the end, matR is stored using a cs column storage  */
//-----------------------------------------------------------------------------------------------------------------------------
void QRSolve(const MatrixEntry& mat, MatrixEntry* matR, VectorEntry* rhs)
{
   if(matR==0 && rhs==0)  return;  //nothing to do !!!
   ValueType vtm=mat.valueType_; StrucType stm=mat.strucType_;
   ValueType vtr=rhs->valueType_;
   //call QRSolve(LargeMatrix,...)
   if(stm == _scalar)
   {
      switch(vtm)
      {
         case _real :
            {
               LargeMatrix<real_t>* lm=0;
               std::vector<real_t>* sm=0;
               if(matR!=0) lm= matR->rEntries_p;
               if(rhs!=0)
               {
                   if(vtr==_real) QRSolve(*mat.rEntries_p, lm, rhs->rEntries_p);
                   else // move matrix to complex because QR does not work yet with real matrix and complex vector
                   {
                       MatrixEntry cmat=mat; cmat.toComplex();
                       LargeMatrix<complex_t>* lmc=0;
                       if(matR!=0)
                         {
                             matR->toComplex();
                             lmc= matR->cEntries_p;
                         }
                       QRSolve(*cmat.cEntries_p, lmc, rhs->cEntries_p);
                       if(matR!=0) matR->toReal(); //come back to real
                   }
               }
               else   QRSolve(*mat.rEntries_p, lm, sm);
            }
            break;
         default : //complex case
            {
               LargeMatrix<complex_t>* lm=0;
               std::vector<complex_t>* sm=0;
               if(matR!=0) lm= matR->cEntries_p;
               if(rhs!=0)
               {
                   if(vtr==_complex) QRSolve(*mat.cEntries_p, lm, rhs->cEntries_p);
                   else
                   {
                       rhs->toComplex();
                       QRSolve(*mat.cEntries_p, lm, rhs->cEntries_p);
                   }
               }
               else   QRSolve(*mat.cEntries_p, lm, sm);
            }
      }
   }
   else
   {
     where("QR(MatrixEntry, ...)");
     error("matrixentry_matrixofmatrices_not_handled");
   }
}

// Gauss solver (only in scalar)
void gaussSolve(MatrixEntry& mat, VectorEntry& b, VectorEntry& x)
{
  ValueType vtm=mat.valueType_;
  ValueType vtb=b.valueType_;
  if(mat.strucType_ != _scalar)
  {
    where("gaussSolve(MatrixEntry, ...)");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(vtm==_real)
  {
    if(vtb==_real)
    {
      x=b;
      gaussSolve(*mat.rEntries_p,*x.rEntries_p);
      return;
    }
    else
    {
      x=b;
      std::vector<std::vector<real_t> > rhss(2);
      rhss[0]=real(*b.cEntries_p);
      rhss[1]=imag(*b.cEntries_p);
      gaussSolve(*mat.rEntries_p,rhss);
      *x.cEntries_p=Vector<real_t>(rhss[0])+complex_t(0,1)*Vector<real_t>(rhss[1]);
      return;
    }
  }
  else
  {
    if(vtb==_real)
    {
      x=b;
      x.toComplex();
      gaussSolve(*mat.cEntries_p,*x.cEntries_p);
      return;
    }
    else
    {
      x=b;
      gaussSolve(*mat.cEntries_p,*x.cEntries_p);
      return;
    }
  }
}

// umfpack solver (only in scalar representation)
// return the reciprocal condition numer
void umfpackSolve(MatrixEntry& mat, VectorEntry& b, VectorEntry& x, real_t& rcond)
{
#ifdef XLIFEPP_WITH_UMFPACK
  ValueType vtm=mat.valueType_;
  ValueType vtb=b.valueType_;
  if(mat.strucType_ != _scalar)
  {
    where("umfpackSolve(MatrixEntry, ...)");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  x=b;  //allocate x with the type of b
  if(vtm==_real)
  {
    if(vtb==_real)
    {
      rcond=mat.rEntries_p->umfpackSolve(*b.rEntries_p, *x.rEntries_p);
      return;
    }
    else
    {
      rcond=mat.rEntries_p->umfpackSolve(*b.cEntries_p, *x.cEntries_p);
      return;
    }
  }
  else
  {
    if(vtb==_real)
    {
      x.toComplex(); // now, x has complex scalar entries
      rcond=mat.cEntries_p->umfpackSolve(*b.rEntries_p, *x.cEntries_p);
      return;
    }
    else
    {
      rcond=mat.cEntries_p->umfpackSolve(*b.cEntries_p, *x.cEntries_p);
      return;
    }
  }
#else
  error("xlifepp_without_umfpack");
#endif // XLIFEPP_WITH_UMFPACK
}

// umfpack solver with multiple right hand side (only in scalar representation)
// xs pointers have to be already allocated
// return the reciprocal condition number
void umfpackSolve(MatrixEntry& mat, std::vector<VectorEntry*>& bs, std::vector<VectorEntry*>& xs, real_t& rcond)
{
#ifdef XLIFEPP_WITH_UMFPACK
  ValueType vtm=mat.valueType_;
  if(mat.strucType_ != _scalar)
  {
    where("umfpackSolve(MatrixEntry, ...)");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  number_t ns=bs.size();
  if(ns==0)
  {
    where("umfpackSolve(MatrixEntry, ...)");
    error("no_rhs");
  }

  //prepare rhs bs
  std::vector<VectorEntry*> cbs=bs;    //copy rhs pointers
  ValueType vtbs=_real;
  for(number_t k=0; k<ns && vtbs==_real; k++) vtbs = bs[k]->valueType_;  // detects if at least one rhs is complex
  if(vtbs==_complex || vtm==_complex) //move real rhs to complex rhs
  {
    for(number_t k=0; k<ns ; k++)
    {
      if(bs[k]->valueType_==_real) //allocate new complex VectorEntry
      {
        cbs[k]= new VectorEntry(complex_t(0), bs[k]->size());
        *cbs[k]+=*bs[k];   //copy real bs[k]
      }
    }
  }
  //prepare results xs
  for(number_t k=0; k<ns; k++) *xs[k]=*cbs[k];

  if(vtm==_real)
  {
    if(vtbs==_real)
    {
      std::vector<std::vector<real_t>* > vb(ns,0), vx(ns,0);
      for(number_t k=0;k<ns; k++)
      {
        vb[k]=bs[k]->rEntries_p;
        vx[k]=xs[k]->rEntries_p;
      }
      rcond=mat.rEntries_p->umfpackSolve(vb,vx);
    }
    else
    {
      std::vector<std::vector<complex_t>* > vb(ns,0), vx(ns,0);
      for(number_t k=0;k<ns; k++)
      {
        vb[k]=bs[k]->cEntries_p;
        vx[k]=xs[k]->cEntries_p;
      }
      rcond=mat.rEntries_p->umfpackSolve(vb,vx);
    }
  }
  else  //complex case
  {
    std::vector<std::vector<complex_t>* > vb(ns,0), vx(ns,0);
    for(number_t k=0;k<ns; k++)
    {
      vb[k]=bs[k]->cEntries_p;
      vx[k]=xs[k]->cEntries_p;
    }
    rcond=mat.cEntries_p->umfpackSolve(vb,vx);
  }

  // deallocate allocated complex rhs
  if(vtbs==_complex || vtm==_complex)
  {
    for(number_t k=0; k<ns ; k++)
    {
      if(bs[k]->valueType_==_real) delete cbs[k]; //deallocate new complex VectorEntry
    }
  }
#else
  error("xlifepp_without_umfpack");
#endif // XLIFEPP_WITH_UMFPACK
}

/*! factorize matrix as LU or LDLt or LDL* along matrix property
   if ft=_noFactorisation, type of factorization is searched
   if matrix is symmetric ldlt factorisation is used
   if matrix is self-adjoint ldlstar factorisation is used
   note that factorisation may be failed because there is no pivoting strategy (except for definite symmetric matrix!)
   if matrix is skew-adjoint ldlstar factorisation may be worked if the diagonal is non zero
   if matrix is skew-symmetric (diagonal is zero!) ldlt failed in any case !
   For the moment LU factorisation is used in case of a skew-adjoint or a skew-symmetric matrix
   in future, for specific algorithm for skew symmetric or skew adjoint matrices see following papers :
    Bunch J.R. Stable Algorithms for Solving Symmetric and Skew-Symmetric Systems, Bull. Austral. Math. Soc., vol 26, 107-119, 1982
    Lau T., Numerical Solution of Skew-Symmetric Linear Systems, Master Thesis, University of British Columbia, 2007

   If umfPack is available, it is used when ft is not specified
   if withPermution = true, LU with row permutation will be used to deal with dense matrices
*/
void factorize(MatrixEntry& A, MatrixEntry& Af, FactorizationType ft, bool withPermutation)
{
  FactorizationType fac = ft;

  if (&Af != &A)
  {
    Af.clear();    //copy A to Af
    Af=A;
  }

  if (ft == _noFactorization) //search factorization type
  {
    fac = _lu;
    if (Af.symmetry() == _symmetric) fac = _ldlt;
    if (Af.symmetry() == _selfAdjoint) fac = _ldlstar;
#ifdef XLIFEPP_WITH_UMFPACK
    fac=_umfpack;
#endif // XLIFEPP_WITH_UMFPACK
  }

  if (fac!=_umfpack && Af.storageType() == _cs) Af.toSkyline(); //convert entries to skyline storage

  //check factorization type and storage consistancy, because we use autofactorization
  if (Af.accessType()==_sym && fac==_lu)
    {//error("structure_not_handled_in_factorization", words("factorization type",_lu), words("access type",_sym));
        fac=_ldlt;
        warning("free_warning","as Matrix has a symmetric storage access, use ldlt instead of lu");
    }

  std::cout<<" factorise matrix "<<Af.nbOfRows()<<" x "<<Af.nbOfCols()<<" using ";
  switch(fac)
  {
    case _ldlt :
    case _llt  :
      std::cout<<"L(D)Lt (skyline) "<<eol;
      Af.ldltFactorize();
      break;
    case _ldlstar :
    case _llstar  :
      std::cout<<"L(D)L* (skyline) "<<eol;
      Af.ldlstarFactorize();
      break;
    case _lu   :
      std::cout<<"LU (skyline) "<<eol;
      Af.luFactorize(withPermutation);
      break;
    case _umfpack   :
      std::cout<<"umfpack (compressed sparse column) "<<eol;
      Af.umfpackFactorize();
      break;
    default :
      error("wrong_factorization_type",words("factorization type",fac));
  }
}

//! factorize matrix as LU or LDLt or LDL*, not preserving A
void factorize(MatrixEntry& A, FactorizationType ft, bool withPermutation) { factorize(A,A,ft,withPermutation); }

void iFactorize(MatrixEntry& A, MatrixEntry& Af, FactorizationType ift)
{
  FactorizationType ifac = ift;

  if (&Af != &A)
  {
    Af.clear();    //copy A to Af
    Af=A;
  }
  if (ift == _noFactorization) //search factorization type
  {
    ifac = _ilu;
    if (Af.symmetry() == _symmetric) ifac = _ildlt;
    if (Af.symmetry() == _selfAdjoint) ifac = _ildlstar;
  }
  std::cout<<" factorise matrix "<<Af.nbOfRows()<<" x "<<Af.nbOfCols()<<" using ";
  switch(ifac)
  {
    case _ildlt :
      std::cout<<"iLDLt (cs) "<<eol;
      Af.ildltFactorize();
      break;
    case _illt  :
      std::cout<<"iLLt (cs) "<<eol;
      Af.illtFactorize();
      break;
    case _ildlstar :
      std::cout<<"iLDL* (cs) "<<eol;
      Af.ildlstarFactorize();
      break;
    case _illstar  :
      std::cout<<"LL* (cs) "<<eol;
      Af.illstarFactorize();
      break;
    case _ilu   :
      std::cout<<"iLU (Cs) "<<eol;
      Af.iluFactorize();
      break;
    default :
      error("wrong_factorization_type",words("factorization type",ifac));
  }
}

//! incomplete factorize matrix as ILU (or iLDLt or iLDL*), not preserving A
void iFactorize(MatrixEntry& A, FactorizationType ift) { iFactorize(A,A,ift); }
//solve linear system when matrix is already factorized
VectorEntry factSolve(MatrixEntry& A, const VectorEntry& B)
{
  trace_p->push("factSolve(MatrixEntry, VectorEntry)");
  if (A.factorization() == _noFactorization) error("term_not_factorized");
  VectorEntry X = B;
  ValueType vtX = _real;
  if(B.valueType_==_complex || A.valueType()==_complex) vtX=_complex;
  if(B.valueType_!=vtX) X.toComplex();
  switch(A.factorization())
  {
    case _ildlt    :
    case _ldlt    :
      A.ldltSolve(B, X);
      break;
    case _ldlstar :
      A.ldlstarSolve(B, X);
      break;
    case _illt      :
    case _llt      :
      A.lltSolve(B, X);
      break;
    case _lu      :
    case _ilu      :
      A.luSolve(B, X);
      break;
    case _umfpack :
      A.umfluSolve(B,X);
      break;
    default :
      error("wrong_factorization_type", words("factorization type",A.factorization()));
  }
  trace_p->pop();
  return X;
}


//-----------------------------------------------------------------------------------------------------------------------------
/*! Resolve an eigen problem with different method
 * This function is a wrapper of different eigenSolvers of class LargeMatrix (e.g: eigenDavidsonSolver, eigenKrylovSchurSolver, ...)
 * \param[in] eigCompMode enumeration of different eigenSolver(davidson, krylovSchur). By default, it should be KrylovSchur
 * \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] sigma shift value used in shift mode. In case of real shift, the real part of sigma will be used
 * \param[in] isRealShift Specify if real shift is used
 * \param[in] isShift Specify if real or complex shift is used */
std::vector<std::pair<complex_t, VectorEntry* > > MatrixEntry::eigenSolve(EigenComputationalMode eigCompMode, number_t nev,
      real_t tol, string_t which, complex_t sigma, bool isRealShift, bool isShift)
{
  real_t realShift = real(sigma);
  complex_t cplxShift = sigma;
  FactorizationType fac = _noFactorization;
  bool isInverted = false;
  LargeMatrix<real_t>* prA = 0;
  LargeMatrix<complex_t>* pcA = 0;
  std::vector<std::pair<complex_t,Vector<complex_t> > > eigenVec;

  if (0 != rEntries_p)  //scalar real matrix
  {
    LargeMatrix<real_t>* pB = 0;
    if (isShift)
    {
      prA = new LargeMatrix<real_t>(_idMatrix, rEntries_p->storageType(), rEntries_p->accessType(),
                                    rEntries_p->nbRows, rEntries_p->nbCols, -realShift);
      prA->toStorage(rEntries_p->storagep());
      *prA += *rEntries_p;
      if (_skyline != prA->storageType())
      {
        LargeMatrix<real_t>* tmp = new LargeMatrix<real_t>(*prA, true);
        tmp->toSkyline();
        delete prA;
        prA = tmp;
      }

      if (_symmetric == prA->sym)
      {
        prA->ldltFactorize();
        fac = _ldlt;
      }
      else
      {
        prA->luFactorize();
        fac = _lu;
      }
      isInverted = true;
    }
    else
    {
      if ("SM" == which)
      {
        prA = new LargeMatrix<real_t>(*rEntries_p, true);
        if (_skyline != prA->storageType()) prA->toSkyline();
        if (_symmetric == prA->sym)
        {
          prA->ldltFactorize();
          fac = _ldlt;
        }
        else
        {
          prA->luFactorize();
          fac = _lu;
        }
        isInverted = true;
      }
      else prA = rEntries_p;
    }

    switch (eigCompMode)
    {
      case _davidson:
        eigenDavidsonSolve<real_t>(prA, pB, eigenVec, nev, tol, which, isInverted, fac, isShift);
        break;
      default:
        eigenKrylovSchurSolve<real_t>(prA, pB, eigenVec, nev, tol, which, isInverted, fac, isShift);
        break;
    }

    if ( prA != rEntries_p) delete prA;
  }

  if (0 != cEntries_p)  //scalar complex matrix
  {
    LargeMatrix<complex_t>* pB = 0;
    if (isShift) //shift != 0
    {
      if (isRealShift) cplxShift = complex_t(realShift,0);
      pcA = new LargeMatrix<complex_t>(_idMatrix, cEntries_p->storageType(), cEntries_p->accessType(),
                                       cEntries_p->nbRows, cEntries_p->nbCols, -cplxShift);
      pcA->toStorage(cEntries_p->storagep());
      *pcA += *cEntries_p;
      if (_skyline != pcA->storageType())
      {
        LargeMatrix<complex_t>* tmp = new LargeMatrix<complex_t>(*pcA, true);
        tmp->toSkyline();
        delete pcA;
        pcA = tmp;
      }

      if (_selfAdjoint == pcA->sym)
      {
        pcA->ldlstarFactorize();
        fac = _ldlstar;
      }
      else if (_symmetric == pcA->sym)
      {
        pcA->ldltFactorize();
        fac = _ldlt;
      }
      else
      {
        pcA->luFactorize();
        fac = _lu;
      }
      isInverted = true;
    }
    else
    {
      if ("SM" == which)
      {
        pcA = new LargeMatrix<complex_t>(*cEntries_p, true);
        if (_skyline != pcA->storageType()) pcA->toSkyline();
        if (_selfAdjoint == pcA->sym)
        {
          pcA->ldlstarFactorize();
          fac = _ldlstar;
        }
        else if (_symmetric == pcA->sym)
        {
          pcA->ldltFactorize();
          fac = _ldlt;
        }
        else
        {
          pcA->luFactorize();
          fac = _lu;
        }
      }
      else pcA = cEntries_p;
    }

    switch (eigCompMode)
    {
      case _davidson :
        eigenDavidsonSolve<complex_t>(pcA, pB, eigenVec, nev, tol, which, isInverted, fac, isShift);
        break;
      default:
        eigenKrylovSchurSolve<complex_t>(pcA, pB, eigenVec, nev, tol, which, isInverted, fac, isShift);
        break;
    }

    if ( pcA != cEntries_p) delete pcA;
  }

  if(0 != rmEntries_p)
  {
    where("eigenSolve(MatrixEntry, ...)");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if(0 != cmEntries_p)
  {
    where("eigenSolve(MatrixEntry, ...)");
    error("matrixentry_matrixofmatrices_not_handled");
  }

  std::vector<std::pair<complex_t, VectorEntry*> > res;
  for (number_t i = 0; i < eigenVec.size(); ++i)
  {
    if (isShift)
    {
      if (0 != rEntries_p) res.push_back(std::make_pair(realShift + 1.0/eigenVec[i].first, new VectorEntry(eigenVec[i].second)));
      else res.push_back(std::make_pair(cplxShift + 1.0/eigenVec[i].first, new VectorEntry(eigenVec[i].second)));
    }
    else
    {
      if ("SM" == which) res.push_back(std::make_pair(1.0/eigenVec[i].first, new VectorEntry(eigenVec[i].second)));
      else res.push_back(std::make_pair(eigenVec[i].first, new VectorEntry(eigenVec[i].second)));
    }
  }

  return (res);
}


//-----------------------------------------------------------------------------------------------------------------------------
/*! Resolve an eigen problem with different method
 * This function is a wrapper of different eigenSolvers of class LargeMatrix (e.g: eigenDavidsonSolver, eigenKrylovSchurSolver, ...)
 * \param[in] pB MatrixEntry containing the matrix B in the eigen problem A*X = lambda*B*X
 * \param[in] eigCompMode Computation mode of eigensolver
 * \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] sigma shift value used in shift mode. In case of real shift, the real part of sigma will be used
 * \param[in] isRealShift Specify if real shift is used
 * \param[in] isShift Specify if real or complex shift is used */
std::vector<std::pair<complex_t, VectorEntry* > > MatrixEntry::eigenSolve(const MatrixEntry* pB, EigenComputationalMode eigCompMode,
      number_t nev, real_t tol, string_t which, complex_t sigma, bool isRealShift, bool isShift)
{
  real_t realShift = real(sigma);
  complex_t cplxShift = sigma;
  FactorizationType fac = _noFactorization;
  bool isInverted = false;
  LargeMatrix<real_t>* prA = rEntries_p, *prB = pB->rEntries_p;
  LargeMatrix<complex_t>* pcA = cEntries_p, *pcB = pB->cEntries_p;
  std::vector<std::pair<complex_t,Vector<complex_t> > > eigenVec;

  if ((0 != rEntries_p) && (0 != pB->rEntries_p))  //scalar real matrix
  {
    if (isShift) // shift != 0
    {
      LargeMatrix<real_t>* tmpB = new LargeMatrix<real_t>(*(pB->rEntries_p));
      *tmpB *= -realShift;
      prA = addMatrixMatrixSkyline(*rEntries_p, *tmpB);

      if (_symmetric == prA->sym)
      {
        prA->ldltFactorize();
        fac = _ldlt;
      }
      else
      {
        prA->luFactorize();
        fac = _lu;
      }
      isInverted = true;
      delete tmpB;
    }
    else
    {
      if ("SM" == which)
      {
        if (_skyline != rEntries_p->storageType()) prA = new LargeMatrix<real_t>(*rEntries_p, true);
        else prA = new LargeMatrix<real_t>(*rEntries_p);
        if (_skyline != prA->storageType()) prA->toSkyline();
        if (_symmetric == prA->sym)
        {
          prA->ldltFactorize();
          fac = _ldlt;
        }
        else
        {
          prA->luFactorize();
          fac = _lu;
        }
        isInverted = true;
      }
      else
      {
        prA = rEntries_p;
        if (_krylovSchur == eigCompMode)
        {
          if (_skyline != pB->rEntries_p->storageType())
          {
            prB = new LargeMatrix<real_t>(*(pB->rEntries_p), true);
            prB->toSkyline();
          }
          else prB = new LargeMatrix<real_t>(*(pB->rEntries_p));

          if (_symmetric == prB->sym)
          {
            prB->ldltFactorize();
            fac = _ldlt;
          }
          else
          {
            prB->luFactorize();
            fac = _lu;
          }
          isInverted = true;
        }
      }
    }

    switch (eigCompMode)
    {
      case _davidson :
        eigenDavidsonSolve<real_t>(prA, prB, eigenVec, nev, tol, which, isInverted, fac, isShift);
        break;
      default:
        eigenKrylovSchurSolve<real_t>(prA, prB, eigenVec, nev, tol, which, isInverted, fac, isShift);
        break;
    }

    if (prA != rEntries_p) delete prA;
    if (prB != pB->rEntries_p) delete prB;
  }

  if ((0 != cEntries_p) && (0 != pB->cEntries_p)) //scalar complex matrix
  {
    if (isShift)
    {
      LargeMatrix<complex_t>* tmpB = new LargeMatrix<complex_t>(*(pB->cEntries_p));
      if (isRealShift) *tmpB *= -realShift;
      else *tmpB *= -cplxShift;
      pcA = addMatrixMatrixSkyline(*cEntries_p, *tmpB);
      if (_selfAdjoint == pcA->sym)
      {
        pcA->ldlstarFactorize();
        fac = _ldlstar;
      }
      else if (_symmetric == pcA->sym)
      {
        pcA->ldltFactorize();
        fac = _ldlt;
      }
      else
      {
        pcA->luFactorize();
        fac = _lu;
      }
      isInverted = true;
      delete tmpB;
    }
    else
    {
      if ("SM" == which)
      {
        if (_skyline != cEntries_p->storageType()) pcA = new LargeMatrix<complex_t>(*cEntries_p, true);
        else pcA = new LargeMatrix<complex_t>(*cEntries_p);
        if (_skyline != pcA->storageType()) pcA->toSkyline();
        if (_selfAdjoint == pcA->sym)
        {
          pcA->ldlstarFactorize();
          fac = _ldlstar;
        }
        else if (_symmetric == pcA->sym)
        {
          pcA->ldltFactorize();
          fac = _ldlt;
        }
        else
        {
          pcA->luFactorize();
          fac = _lu;
        }
        isInverted = true;
      }
      else
      {
        pcA = cEntries_p;
        if (_krylovSchur == eigCompMode)
        {
          if (_skyline != pB->cEntries_p->storageType())
          {
            pcB = new LargeMatrix<complex_t>(*(pB->cEntries_p), true);
            pcB->toSkyline();
          }
          else pcB = new LargeMatrix<complex_t>(*(pB->cEntries_p));

          if (_selfAdjoint == pcB->sym)
          {
            pcB->ldlstarFactorize();
            fac = _ldlstar;
          }
          else if (_symmetric == pcB->sym)
          {
            pcB->ldltFactorize();
            fac = _ldlt;
          }
          else
          {
            pcB->luFactorize();
            fac = _lu;
          }
          isInverted = true;
        }
      }
    }

    switch (eigCompMode)
    {
      case _davidson :
        eigenDavidsonSolve<complex_t>(pcA, pcB, eigenVec, nev, tol, which, isInverted, fac, isShift);
        break;
      default:
        eigenKrylovSchurSolve<complex_t>(pcA, pcB, eigenVec, nev, tol, which, isInverted, fac, isShift);
        break;
    }

    if (pcA != cEntries_p) delete pcA;
    if (pcB != pB->cEntries_p) delete pcB;
  }

  if((0 != rmEntries_p) || (0 != pB->rmEntries_p))
  {
    where("eigenSolve(MatrixEntry, ...)");
    error("matrixentry_matrixofmatrices_not_handled");
  }
  if((0 != cmEntries_p) || (0 != pB->cmEntries_p))
  {
    where("eigenSolve(MatrixEntry, ...)");
    error("matrixentry_matrixofmatrices_not_handled");
  }

  std::vector<std::pair<complex_t, VectorEntry*> > res;
  for (number_t i = 0; i < eigenVec.size(); ++i)
  {
     if (isShift)
    {
      if (0 != rEntries_p) res.push_back(std::make_pair(realShift + 1.0/eigenVec[i].first, new VectorEntry(eigenVec[i].second)));
      else res.push_back(std::make_pair(cplxShift + 1.0/eigenVec[i].first,new VectorEntry(eigenVec[i].second)));
    }
    else
    {
      if ("SM" == which) res.push_back(std::make_pair(1.0/eigenVec[i].first, new VectorEntry(eigenVec[i].second)));
      else res.push_back(std::make_pair(eigenVec[i].first, new VectorEntry(eigenVec[i].second)));
    }
  }

  return (res);
}

} //end of namespace xlifepp