xref: /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/largeMatrix/MatrixStorage.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 MatrixStorage.cpp
  \authors D. Martin, E. Lunéville, M.-H. Nguyen, N. Kielbasiewicz
  \since 21 jun 2011
  \date 21 feb 2013

  \brief Implementation of xlifepp::MatrixStorage class and functionnalities
*/

#include "MatrixStorage.hpp"
#include "denseStorage/ColDenseStorage.hpp"
#include "denseStorage/RowDenseStorage.hpp"
#include "denseStorage/DualDenseStorage.hpp"
#include "denseStorage/SymDenseStorage.hpp"
#include "csStorage/ColCsStorage.hpp"
#include "csStorage/RowCsStorage.hpp"
#include "csStorage/DualCsStorage.hpp"
#include "csStorage/SymCsStorage.hpp"
#include "skylineStorage/SkylineStorage.hpp"
#include "skylineStorage/DualSkylineStorage.hpp"
#include "skylineStorage/SymSkylineStorage.hpp"

namespace xlifepp
{
/*---------------------------------------------------------------------------------------
  Constructors, destructor
  ---------------------------------------------------------------------------------------*/
// default constructor
MatrixStorage::MatrixStorage(string_t id)
  : storageType_(_noStorage), accessType_(_noAccess), nbRows_(0), nbCols_(0), nbObjectsSharingThis_(0), stringId(id) {}

// construct MatrixStorage with no size
MatrixStorage::MatrixStorage(StorageType st, AccessType at, string_t id)
  : storageType_(st), accessType_(at), nbRows_(0), nbCols_(0), nbObjectsSharingThis_(0), stringId(id)
{
  theMatrixStorages.push_back(this);
}

// construct MatrixStorage with size
MatrixStorage::MatrixStorage(StorageType st, AccessType at, number_t nr, number_t nc, string_t id)
  : storageType_(st), accessType_(at), nbRows_(nr), nbCols_(nc), nbObjectsSharingThis_(0), stringId(id)
{
  if(trackingObjects) theMatrixStorages.push_back(this);
}

//copy constructor, NOTE : nbObjectsSharingThis_ is reset to 0
MatrixStorage::MatrixStorage(const MatrixStorage& ms)
{
  storageType_=ms.storageType_;
  accessType_=ms.accessType_;
  nbRows_=ms.nbRows_;
  nbCols_=ms.nbCols_;
  nbObjectsSharingThis_=0;
}

// create matrix storage from indices (vector of column indices for each row)
// like a constructor
// if storage already exists we force its creation by changing id
MatrixStorage* createMatrixStorage(StorageType st, AccessType at, number_t nbr, number_t nbc,
                                   const std::vector<std::vector<number_t> >& indices, const string_t & idu)
{
  string_t id=idu;
  number_t k=0;
  MatrixStorage* ms= findMatrixStorage(id,st,at);
  while(ms!=0)
    {
      id=idu+"_"+tostring(k);
      ms= findMatrixStorage(id,st,at);
      k++;
    }
  switch(st)
    {
      case _cs :
        switch(at)
          {
            case _row  : return new RowCsStorage(nbr, nbc, indices, id);
            case _col  :
              {
                std::vector<std::vector<number_t> > rowindices(nbc);
                std::vector<std::vector<number_t> >::const_iterator itr=indices.begin();
                std::vector<number_t>::const_iterator itc;
                number_t r=1;
                for(; itr!=indices.end(); itr++, r++)
                  for(itc=itr->begin(); itc!= itr->end(); itc++)
                    rowindices[*itc - 1].push_back(r);
                return new ColCsStorage(nbr, nbc, rowindices, id);
              }
            case _dual : return new DualCsStorage(nbr, nbc, indices, id);
            case _sym  : return new SymCsStorage(nbr, indices,_all, id);   //nbr=nbc if symmetric
            default :
              where("createMatrixStorage");
              error("storage_bad_access",words("access type",at),words("storage type",st));
          }
        break;
      case _skyline :
        switch(at)
          {
            case _dual : return new DualSkylineStorage(nbr, nbc, indices, id);
            case _sym  : return new SymSkylineStorage(nbr, indices, id); //nbr=nbc if symmetric
            default :
              where("createMatrixStorage");
              error("storage_bad_access",words("access type",at),words("storage type",st));
          }
      case _dense :
        switch(at)
          {
            case _row  : return new RowDenseStorage(nbr, nbc, id);
            case _col  : return new ColDenseStorage(nbr, nbc, id);
            case _dual : return new DualDenseStorage(nbr, nbc, id);
            case _sym  : return new SymDenseStorage(nbr, id); //nbr=nbc if symmetric
            default :
              where("createMatrixStorage");
              error("storage_bad_access",words("access type",at),words("storage type",st));
          }
      default :
        where("createMatrixStorage");
        error("storage_not_handled",words("storage type",st),words("storage type",at));
    }
  return 0;
}

// construct a storage from submatrix of current storage, storage type is not changed
MatrixStorage* MatrixStorage::extract(const std::vector<number_t>& rowIndex, const std::vector<number_t>& colIndex, const string_t& ids)
{
  number_t nbr=rowIndex.size(), nbc = colIndex.size();
  if (nbr==0 || nbr>nbRows_)
  {
    where("MatrixStorage::extract(...)");
    error("dim_not_in_range", 1, nbRows_);
  }
  if (nbc==0 || nbc>nbCols_)
  {
    where("MatrixStorage::extract(...)");
    error("dim_not_in_range", 1, nbCols_);
  }
  string_t id= ids;
  if(id=="") id=stringId+"_submatrix";

  if(storageType_==_dense) //dense storage
    {
      switch(accessType_)
        {
          case _row  : return new RowDenseStorage(nbr, nbc, id);
          case _col  : return new ColDenseStorage(nbr, nbc, id);
          case _dual : return new DualDenseStorage(nbr, nbc, id);
          case _sym  : return new SymDenseStorage(nbr, id); //nbr=nbc if symmetric
          default :
            where("MatrixStorage constructor");
            error("storage_bad_access",words("access type",accessType_),words("storage type",storageType_));
        }
    }

  //extract indices array
  std::vector<number_t> revCol(nbCols_,0);  //reverse colIndex
  std::vector<number_t>::const_iterator itn=colIndex.begin();
  for(number_t i=1; i<=nbc; ++i, ++itn) revCol[*itn]=i;
  std::vector<std::vector<number_t> > indices(nbr);
  itn=rowIndex.begin();
  for(number_t r=0; r<nbr; r++, ++itn)
    {
      std::vector<std::pair<number_t, number_t> > rowr=getRow(_noSymmetry, *itn);
      std::vector<std::pair<number_t, number_t> > ::iterator itp=rowr.begin();
      number_t s=0;
      for(; itp!=rowr.end(); ++itp)
        if(revCol[itp->first]>0) s++;
      if(s>0)
        {
          std::vector<number_t>& indr=indices[r];
          indr.reserve(s);
          for(; itp!=rowr.end(); ++itp)
            {
              number_t c=revCol[itp->first];
              if(c>0) indr.push_back(c);
            }
        }
    }
  //general construction for cs or skyline, could be improved for skyline
  return createMatrixStorage(storageType_, accessType_, nbr, nbc, indices, id);
}

// destruct MatrixStorage, error if storage is shared by other matrices
MatrixStorage::~MatrixStorage()
{
  if(nbObjectsSharingThis_ > 0)
    {
      error("storage_baddelete", name(), nbObjectsSharingThis_);
    }
  // remove pointer of current object from vector of MatrixStorage*
  std::vector<MatrixStorage*>::iterator it(find(theMatrixStorages.begin(), theMatrixStorages.end(), this));
  if(it != theMatrixStorages.end())
    {
      theMatrixStorages.erase(it);
    }
}

// delete all MatrixStorage objects
void MatrixStorage::clearGlobalVector()
{
  while(MatrixStorage::theMatrixStorages.size() > 0)
    {
      MatrixStorage* ms = MatrixStorage::theMatrixStorages[0];
      if(ms!=0)
        {
          if(ms->nbObjectsSharingThis_>0)  //force deletion
            {
              warning("storage_hasardousdelete", ms->name(), ms->nbObjectsSharingThis_);
              ms->nbObjectsSharingThis_=0;
            }
          delete MatrixStorage::theMatrixStorages[0];
        }
    }
}

//--------------------------------------------------------------------------
//  Error handlers
//--------------------------------------------------------------------------
void MatrixStorage::noFactorization(const string_t& f) const
{
  theMessageData << f + " " + name() + " no factorization";
  error("largematrix_nofactorization", theMessageData);
}

void MatrixStorage::noSolver(const string_t& f) const
{
  theMessageData << f + " " + name() + " no solver";
  error("largematrix_nosolver", theMessageData);
}

void MatrixStorage::isSingular(const string_t& f, const number_t r) const
{
  theMessageData << f + " " + name() + " is singular" << r;
  error("largematrix_singular", theMessageData);
}

/*--------------------------------------------------------------------------------
   "visualization" to output stream of the structure of matrix storage
   using virtual function pos(i,j), pos(i,j) =0 means that (i,j) is outside the storage
   part must be equal to _all, _lower or _upper
--------------------------------------------------------------------------------*/
void MatrixStorage::visual(std::ostream& os) const
{
  if(theVerboseLevel > 0)
    {
      printHeader(os);
    }
  if(theVerboseLevel < 2)
    {
      return;
    }
  number_t rmax = std::min(nbRows_, number_t(10 * theVerboseLevel)),
           cmax = std::min(nbCols_, number_t(10 * theVerboseLevel + 5));
  os << std::setw(11);
  for(number_t j = 1; j <= cmax; j++)
    {
      os << j % 10;    //col numbers
    }
  for(number_t r = 1; r <= rmax; r++)
    {
      string_t str(cmax, '.');
      for(number_t c = 1; c <= cmax; c++)
        {
          number_t p = pos(r, c);
          if(p != 0 && (accessType_ != _sym || (accessType_ == _sym && r > c)))
            {
              str.at(c - 1) = 'x';
            }
          if(r == c)
            {
              str.at(c - 1) = 'd';
            }
        }
      os << std::endl << std::setw(8) << r << "  " << str;
      if(cmax < nbCols_)
        {
          os << " ...(continued)";
        }
    }
  os << std::endl << std::setw(11);
  for(number_t j = 1; j <= cmax; j++)
    {
      os << j % 10;
    }
  os << std::endl;
}

/*--------------------------------------------------------------------------------
   I/O utilities
--------------------------------------------------------------------------------*/
void MatrixStorage::printHeader(std::ostream& os) const
{
  os << words("matrix storage") << " :" << name() << ", " << words("size") << " = " << size()
     << ", id = " << stringId
     << ", " << nbRows_ << " " << words("rows") << ", " << nbCols_ << " " << words("columns")
     << " (" << words("shared by") << " " << nbObjectsSharingThis_ << " " << words("objects") << ").\n";
}

void MatrixStorage::print(std::ostream& os) const
{
  printHeader(os);
  //see child if something to write
}

std::ostream& operator<<(std::ostream& os, const MatrixStorage& ms)
{
  ms.print(os);
  return os;
}

void printCoo(std::ostream& os, const real_t& v, number_t i, number_t j, real_t tol)
{
  if(std::abs(v) > tol)
    {
      os << i << " " << j << " " << v << std::endl;
    }
}

void printCoo(std::ostream& os, const complex_t& v, number_t i, number_t j, real_t tol)
{
  if(std::abs(v) > tol)
    {
      os << i << " " << j << " " << v.real() << " " << v.imag() << std::endl;
    }
}

void printDense(std::ostream& os, const real_t& v, number_t k)
{
  os << v << " ";
}

void printDense(std::ostream& os, const complex_t& v, number_t k)
{
  os << v.real() << " " << v.imag() << "   ";
}

void readItem(std::istream& ifs, complex_t& v, bool realAsCmplx)
{
  real_t r1 = 0.;
  real_t r2 = 0.;
  if(realAsCmplx)
    {
      ifs >> r1;    //cast real to complex
    }
  else
    {
      ifs >> r1 >> r2;
    }
  v = complex_t(r1, r2);
}

void readItem(std::istream& ifs, real_t& v, bool realAsCmplx)
{
  ifs >> v;
}

number_t numberOfRows(const real_t& v)
{
  return 1;
}
number_t numberOfRows(const complex_t& v)
{
  return 1;
}
number_t numberOfCols(const real_t& v)
{
  return 1;
}
number_t numberOfCols(const complex_t& v)
{
  return 1;
}

// find matrix storage in vector theMatrixStorages of class MatrixStorage
MatrixStorage* findMatrixStorage(const string_t& id, StorageType st, AccessType at)
{
  for(std::vector<MatrixStorage*>::iterator it = MatrixStorage::theMatrixStorages.begin(); it != MatrixStorage::theMatrixStorages.end(); it++)
    {
      if(id == (*it)->stringId && (*it)->storageType()==st && (*it)->accessType()==at)
        {
          return *it;
        }
    }
  return 0;
}

void printMatrixStorages(std::ostream& out)
{
  for(std::vector<MatrixStorage*>::iterator it = MatrixStorage::theMatrixStorages.begin(); it != MatrixStorage::theMatrixStorages.end(); it++)
    {
      (*it)->printHeader(out);
    }
}

//! get (row indices, adress) of col c in set [r1,r2], generic algorithm, see child for better algorithms
std::vector<std::pair<number_t, number_t> > MatrixStorage::getCol(SymType s, number_t c, number_t r1, number_t r2) const
{
  if(accessType_!=_sym) s=_noSymmetry;
  number_t lr=r2;
  std::vector<std::pair<number_t, number_t> > rowadrs;
  if(lr==0) lr=nbRows_;
  if(lr<r1) return rowadrs;
  number_t nbr=lr-r1+1;
  rowadrs.resize(nbr);
  std::vector<number_t> cols(1,c), rows(nbr,1);
  std::vector<number_t>::iterator itr=rows.begin();
  for(number_t r=r1; r<=lr; r++,itr++) *itr=r;
  std::vector<number_t> adrs;
  positions(rows,cols,adrs,false,s);
  std::vector<std::pair<number_t, number_t> >::iterator itra=rowadrs.begin();
  itr=adrs.begin();
  number_t k=0;
  for(number_t r=r1; r<=lr; r++,itr++)
    if(*itr!=0)
      {
        *itra=std::make_pair(r,*itr);
        itra++; k++;
      }
  rowadrs.resize(k);
  return rowadrs;
}

std::vector<std::pair<number_t, number_t> > MatrixStorage::getRow(SymType s, number_t r, number_t c1, number_t c2) const
{
  if(accessType_!=_sym) s=_noSymmetry;
  number_t lc=c2;
  std::vector<std::pair<number_t, number_t> > coladrs;
  if(lc==0) lc=nbCols_;
  if(lc<c1) return coladrs;
  number_t nbc=lc-c1+1;
  coladrs.resize(nbc);
  std::vector<number_t> rows(1,r), cols(nbc,1);
  std::vector<number_t>::iterator itc=cols.begin();
  for(number_t c=c1; c<=lc; c++,itc++) *itc=c;
  std::vector<number_t> adrs;
  positions(rows,cols,adrs,false,s);
  std::vector<std::pair<number_t, number_t> >::iterator itca=coladrs.begin();
  itc=adrs.begin();
  number_t k=0;
  for(number_t c=c1; c<=lc; c++,itc++,itca++)
    if(*itc!=0)
      {
        *itca=std::make_pair(c,*itc);
        k++;
      }
  coladrs.resize(k);
  return coladrs;
}


//! get col indices of row r in set [c1,c2], generic algorithm using position, see children for better algorithms
std::set<number_t> MatrixStorage::getCols(number_t r, number_t c1, number_t c2) const
{
  number_t lc=c2;
  std::set<number_t> colset;
  if(lc==0) lc=nbCols_;
  if(lc<c1) return colset;
  number_t nbc=lc-c1+1;
  std::vector<number_t> rows(1,r), cols(nbc,1);
  std::vector<number_t>::iterator itc=cols.begin();
  for(number_t c=c1; c<=lc; c++,itc++) *itc=c;
  std::vector<number_t> adrs;
  positions(rows,cols,adrs,false);  //retry adress
  itc=adrs.begin();
  for(number_t c=c1; c<=lc; c++, itc++)
    if(*itc!=0) colset.insert(c);
  return colset;
}

//! get row indices of col c in set [r1,r2], generic algorithm using position, see children for better algorithms
std::set<number_t> MatrixStorage::getRows(number_t c, number_t r1, number_t r2) const
{
  number_t lr=r2;
  std::set<number_t> rowset;
  if(lr==0) lr=nbRows_;
  if(lr<r1) return rowset;
  number_t nbr=lr-r1+1;
  std::vector<number_t> cols(1,c), rows(nbr,1);
  std::vector<number_t>::iterator itr=rows.begin();
  for(number_t r=r1; r<=lr; r++,itr++) *itr=r;
  std::vector<number_t> adrs;
  positions(rows,cols,adrs,false);  //retry adress
  itr=adrs.begin();
  for(number_t r=r1; r<=lr; r++, itr++)
    if(*itr!=0) rowset.insert(r);
  return rowset;
}

//! get diagonal adresses, algorithm using position in non sym/dual storage cases
std::vector<number_t> MatrixStorage::getDiag() const
{
  if(accessType_==_sym || accessType_==_dual)       //diagonal is stored at the beginning
    return trivialNumbering<number_t>(1,diagonalSize());
  //other cases
  number_t nd=diagonalSize();
  std::vector<number_t> diag(nd);
  std::vector<number_t>::iterator itd=diag.begin();
  for(number_t i=1; i<=nd; i++, itd++) *itd=pos(i,i);
  return diag;
}

/*!
  get col indices of row r in set [c1,c2], generic algorithm using position, see children for better algorithms
  update a vector of numbers with no reallocation, size it before. nbcol is the number of loaded col indices
*/
void MatrixStorage::getColsV(std::vector<number_t>& colV, number_t& nbcol, number_t r, number_t c1, number_t c2) const
{
  number_t lc=c2;
  if(lc==0) lc=nbCols_;
  if(lc<c1) {nbcol=0; return;};
  number_t nbc=lc-c1+1;
  std::vector<number_t> rows(1,r), cols(nbc,1);
  std::vector<number_t>::iterator itc=cols.begin();
  for(number_t c=c1; c<=lc; c++,itc++) *itc=c;
  std::vector<number_t> adrs;
  positions(rows,cols,adrs,false);  //retry adress
  itc=adrs.begin();
  std::vector<number_t>::iterator itV=colV.begin();
  for(number_t c=c1; c<=lc; c++, itc++)
    if(*itc!=0) {*itV=c; ++itV; nbcol++;}
}

/*!
  get col indices of row r in set [c1,c2], generic algorithm using position, see children for better algorithms
  update a vector of numbers with no reallocation, size it before. nbrow is the number of loaded row indices
*/
void MatrixStorage::getRowsV(std::vector<number_t>& rowV, number_t& nbrow, number_t c, number_t r1, number_t r2) const
{
  number_t lr=r2;
  if(lr==0) lr=nbCols_;
  nbrow=0;
  if(lr<r1) return;
  number_t nbr=lr-r1+1;
  std::vector<number_t> rows(1,nbr), cols(1,c);
  std::vector<number_t>::iterator itr=rows.begin();
  for(number_t r=r1; r<=lr; r++,itr++) *itr=r;
  std::vector<number_t> adrs;
  positions(rows,cols,adrs,false);  //retry adress
  itr=adrs.begin();
  std::vector<number_t>::iterator itV=rowV.begin();
  for(number_t r=r1; r<=lr; r++, itr++)
    if(*itr!=0) {*itV=r; ++itV; nbrow++;}
}

/*! add storage in current storage with row/col mapping (generic algorithm, use addRow, addCol, getRow and getCol of children)
    ms is the matrix storage to add to current storage
    rowmap(r) gives the row number in current storage of row r of MatrixStorage ms (start index : 1)
    colmap(c) gives the col number in current storage of col c of MatrixStorage ms (start index : 1)
    if rowmap (resp. colmap) is void, it means that row numberings are the same for the both storages
    !!! warning : ms storage has to be smaller than the current storage
*/
void MatrixStorage::add(const MatrixStorage& ms, const std::vector<number_t>& rowmap, const std::vector<number_t>& colmap)
{
  if(storageType_==_dense) return;   //nothing to do

  std::vector<number_t> rowmap_=rowmap, colmap_=colmap;
  std::vector<number_t>::iterator it;
  if(rowmap_.size()==0)
    {
      rowmap_.resize(nbRows_);
      number_t k=1;
      for(it=rowmap_.begin(); it!=rowmap_.end(); it++, k++) *it=k;
    }
  if(colmap_.size()==0)
    {
      colmap_.resize(nbCols_);
      number_t k=1;
      for(it=colmap_.begin(); it!=colmap_.end(); it++, k++) *it=k;
    }
  it=std::max_element(rowmap_.begin(),rowmap_.end());
  if(*it>nbRows_)
    {
      where("MatrixStorage::add");
      error("storage_numbering_too_large",words("row"),"ms");
    }
  it=std::max_element(colmap_.begin(),colmap_.end());
  if(*it>nbCols_)
    {
      where("MatrixStorage::add");
      error("storage_numbering_too_large",words("column"),"ms");
    }

  if(accessType_==_row || accessType_==_dual || accessType_==_sym) //row algorithm
    {
      MatrixPart mp=_all;
      if(accessType_!=_row) mp=_lower;
      number_t k=1, n=ms.nbCols_;
      for(it=rowmap_.begin(); it!=rowmap_.end(); it++, k++)
        {
          //if(accessType_==_dual) n=k;
          std::set<number_t> cols=ms.getCols(k,1,n);          // retry col indices of row k in ms storage
          std::set<number_t>::iterator its=cols.begin();
          std::set<number_t> mapcols;
          for(; its!=cols.end(); its++) mapcols.insert(colmap_[*its-1]); // apply col numbering map
          addRow(*it,mapcols,mp);
        }
    }
  if(accessType_==_col || accessType_==_dual) //col algorithm
    {
      MatrixPart mp=_all;
      if(accessType_!=_col) mp=_upper;
      number_t k=1,  n=ms.nbRows_;
      for(it=colmap_.begin(); it!=colmap_.end(); it++, k++)
        {
          //if(accessType_==_dual) n=k;
          std::set<number_t> rows=ms.getRows(k,1,n);           // retry row indices of col k in ms storage
          std::set<number_t>::iterator its=rows.begin();
          std::set<number_t> maprows;
          for(; its!=rows.end(); its++) maprows.insert(rowmap_[*its-1]);  // apply col numbering map
          addCol(*it,maprows,mp);
        }
    }
  return;
}

/*! create scalar col indices for each row from current storage and submatrix sizes
    this is a generic tool used by toScalar() member function of storage children
    it is based on getRow function (virtual function)
    nbr : row size of submatrix
    nbc : col size of submatrix
    if nbc=nbr=1 (scalar case) return the same col indices of current storage
        See children for optimized algorithms
*/
std::vector<std::vector<number_t> > MatrixStorage::scalarColIndices(dimen_t nbr, dimen_t nbc)
{
  std::vector<std::vector<number_t> > cols(nbr*nbRows_);
  std::vector<std::vector<number_t> >::iterator itcs=cols.begin();
  for(number_t r=0; r<nbRows_; r++)
    {
      std::vector<std::pair<number_t, number_t> > row = getRow(_noSymmetry,r+1);
      std::vector<std::pair<number_t, number_t> >::iterator itr;
      for(dimen_t rr=0; rr<nbr; rr++, itcs++)
        {
          itcs->resize(nbc*row.size());
          std::vector<number_t>::iterator itc=itcs->begin();
          for(itr=row.begin(); itr!=row.end(); itr++)
            for(dimen_t cc=0; cc <nbc; cc++, itc++) *itc= nbc*(itr->first-1) + cc + 1;
        }
    }
  return cols;
}

// Matrix * Matrix (result in row dense storage) - generic algorithm - (see child bor better algorithm)
void MatrixStorage::multMatrixMatrix(const std::vector<real_t>& mA, const MatrixStorage& stB, const std::vector<real_t>& mB, std::vector<real_t>& mR, SymType symA, SymType symB) const
{
  multMatrixMatrixGeneric(mA.begin(), stB, mB.begin(), mR.begin(),symA, symB);
}
void MatrixStorage::multMatrixMatrix(const std::vector<real_t>& mA, const MatrixStorage& stB, const std::vector<complex_t>& mB, std::vector<complex_t>& mR, SymType symA, SymType symB) const
{
  multMatrixMatrixGeneric(mA.begin(), stB, mB.begin(), mR.begin(),symA,symB);
}
void MatrixStorage::multMatrixMatrix(const std::vector<complex_t>& mA, const MatrixStorage& stB, const std::vector<real_t>& mB, std::vector<complex_t>& mR, SymType symA, SymType symB) const
{
  multMatrixMatrixGeneric(mA.begin(), stB, mB.begin(), mR.begin(),symA, symB);
}
void MatrixStorage::multMatrixMatrix(const std::vector<complex_t>& mA, const MatrixStorage& stB, const std::vector<complex_t>& mB, std::vector<complex_t>& mR, SymType symA, SymType symB) const
{
  multMatrixMatrixGeneric(mA.begin(), stB, mB.begin(), mR.begin(),symA, symB);
}

//matrix of matrices
void MatrixStorage::multMatrixMatrix(const std::vector<Matrix<real_t> >& mA, const MatrixStorage& stB, const std::vector<Matrix<real_t> >& mB, std::vector<Matrix<real_t> >& mR, SymType symA, SymType symB) const
{
  multMatrixMatrixGeneric(mA.begin(), stB, mB.begin(), mR.begin(),symA, symB);
}
void MatrixStorage::multMatrixMatrix(const std::vector<Matrix<real_t> >& mA, const MatrixStorage& stB, const std::vector<Matrix<complex_t> >& mB, std::vector<Matrix<complex_t> >& mR, SymType symA, SymType symB) const
{
  multMatrixMatrixGeneric(mA.begin(), stB, mB.begin(), mR.begin(),symA, symB);
}
void MatrixStorage::multMatrixMatrix(const std::vector<Matrix<complex_t> >& mA, const MatrixStorage& stB, const std::vector<Matrix<real_t> >& mB, std::vector<Matrix<complex_t> >& mR, SymType symA, SymType symB) const
{
  multMatrixMatrixGeneric(mA.begin(), stB, mB.begin(), mR.begin(),symA, symB);
}
void MatrixStorage::multMatrixMatrix(const std::vector<Matrix<complex_t> >& mA, const MatrixStorage& stB, const std::vector<Matrix<complex_t> >& mB, std::vector<Matrix<complex_t> >& mR, SymType symA, SymType symB) const
{
  multMatrixMatrixGeneric(mA.begin(), stB, mB.begin(), mR.begin(),symA, symB);
}

// Matrix * diag Matrix (result in current storage) - generic algorithm - (see child bor better algorithm)
void MatrixStorage::multMatrixDiagMatrix(const std::vector<real_t>& mA, const std::vector<real_t>& vB, std::vector<real_t>& mR) const
{
  multMatrixDiagMatrixGeneric(mA.begin(), vB.begin(), mR.begin());
}
void MatrixStorage::multMatrixDiagMatrix(const std::vector<real_t>& mA, const std::vector<complex_t>& vB, std::vector<complex_t>& mR) const
{
  multMatrixDiagMatrixGeneric(mA.begin(), vB.begin(), mR.begin());
}
void MatrixStorage::multMatrixDiagMatrix(const std::vector<complex_t>& mA, const std::vector<real_t>& vB, std::vector<complex_t>& mR) const
{
  multMatrixDiagMatrixGeneric(mA.begin(), vB.begin(), mR.begin());
}
void MatrixStorage::multMatrixDiagMatrix(const std::vector<complex_t>& mA, const std::vector<complex_t>& vB, std::vector<complex_t>& mR) const
{
  multMatrixDiagMatrixGeneric(mA.begin(), vB.begin(), mR.begin());
}
//matrix of matrices
void MatrixStorage::multMatrixDiagMatrix(const std::vector<Matrix<real_t> >& mA, const std::vector<Vector<real_t> >& vB, std::vector<Matrix<real_t> >& mR)const
{
  multMatrixDiagMatrixGeneric(mA.begin(), vB.begin(), mR.begin());
}
void MatrixStorage::multMatrixDiagMatrix(const std::vector<Matrix<real_t> >& mA, const std::vector<Vector<complex_t> >& vB, std::vector<Matrix<complex_t> >& mR)const
{
  multMatrixDiagMatrixGeneric(mA.begin(), vB.begin(), mR.begin());
}
void MatrixStorage::multMatrixDiagMatrix(const std::vector<Matrix<complex_t> >& mA, const std::vector<Vector<real_t> >& vB, std::vector<Matrix<complex_t> >& mR) const
{
  multMatrixDiagMatrixGeneric(mA.begin(), vB.begin(), mR.begin());
}
void MatrixStorage::multMatrixDiagMatrix(const std::vector<Matrix<complex_t> >& mA, const std::vector<Vector<complex_t> >& vB, std::vector<Matrix<complex_t> >& mR)const
{
  multMatrixDiagMatrixGeneric(mA.begin(), vB.begin(), mR.begin());
}

// diag Matrix * Matrix (result in current storage) - generic algorithm - (see child bor better algorithm)
void MatrixStorage::multDiagMatrixMatrix(const std::vector<real_t>& vA, const std::vector<real_t>& mB, std::vector<real_t>& mR) const
{
  multDiagMatrixMatrixGeneric(vA.begin(), mB.begin(), mR.begin());
}
void MatrixStorage::multDiagMatrixMatrix(const std::vector<real_t>& vA, const std::vector<complex_t>& mB, std::vector<complex_t>& mR) const
{
  multDiagMatrixMatrixGeneric(vA.begin(), mB.begin(), mR.begin());
}
void MatrixStorage::multDiagMatrixMatrix(const std::vector<complex_t>& vA, const std::vector<real_t>& mB, std::vector<complex_t>& mR) const
{
  multDiagMatrixMatrixGeneric(vA.begin(), mB.begin(), mR.begin());
}
void MatrixStorage::multDiagMatrixMatrix(const std::vector<complex_t>& vA, const std::vector<complex_t>& mB, std::vector<complex_t>& mR) const
{
  multDiagMatrixMatrixGeneric(vA.begin(), mB.begin(), mR.begin());
}

void MatrixStorage::multDiagMatrixMatrix(const std::vector<Vector<real_t> >& vA, const std::vector<Matrix<real_t> >& mB, std::vector<Matrix<real_t> >& mR)const
{
  multDiagMatrixMatrixGeneric(vA.begin(), mB.begin(), mR.begin());
}
void MatrixStorage::multDiagMatrixMatrix(const std::vector<Vector<real_t> >& vA, const std::vector<Matrix<complex_t> >& mB, std::vector<Matrix<complex_t> >& mR)const
{
  multDiagMatrixMatrixGeneric(vA.begin(), mB.begin(), mR.begin());
}
void MatrixStorage::multDiagMatrixMatrix(const std::vector<Vector<complex_t> >& vA, const std::vector<Matrix<real_t> >& mB, std::vector<Matrix<complex_t> >& mR) const
{
  multDiagMatrixMatrixGeneric(vA.begin(), mB.begin(), mR.begin());
}
void MatrixStorage::multDiagMatrixMatrix(const std::vector<Vector<complex_t> >& vA, const std::vector<Matrix<complex_t> >& mB, std::vector<Matrix<complex_t> >& mR)const
{
  multDiagMatrixMatrixGeneric(vA.begin(), mB.begin(), mR.begin());
}

/* =====================================================================================
   colamd/symamd - sparse matrix column ordering algorithm adapt from SuiteSparse colamd
  ======================================================================================
  from compressed sparse storage data (rowIndices, colPointer) produces
  column permutation Q (colPermvector) of A such that P(AQ)=LU or (AQ)'AQ=LL' have less fill-in
  return false if fails

  nbr        : number of rows
  nbc        : number of columns
  rowIndices : list of row indices of non zeros
  colPointer : pointer to column entries
  colPerm    : the vector permutation of columns (size nbc)

     exemple : nbr=5, nbc=4 matrix with 11 nonzero entries
            row/col  0 1 2 3
               0     x 0 x 0
               1     x 0 x x
               2     0 x x 0
               3     0 0 x x
               4     x x 0 0
            rowIndices = {0, 1, 4,  2, 4,  0, 1, 2, 3,  1, 3} ;
            colPinter  = {0,        3,     5,           9,   11} ;
    note : nnz = colPointer[nbc]= rowIndices.size()

*/
bool colAmd(number_t nbr, number_t nbc, const std::vector<number_t>& rowIndices,
            const std::vector<number_t>& colPointer, std::vector<number_t>& colPerm)
{
//
//   // copy rowIndices and resize to  2*nnz + n_col + Col_size + Row_size in rowIndices
//   number_t nnz=rowIndices.size(), s= 2*nnz + n_col + Col_size + Row_size
//   std::vector<number_t> A(rowIndices)
  error("not_yet_implemented","colAmd(...)");
  return false; // dummy return
}

bool symAmd(number_t nbr, number_t nbc, const std::vector<number_t>& colIndices,
            const std::vector<number_t>& rowPointer, std::vector<number_t>& colPerm)
{
  error("not_yet_implemented","symAmd(...)");
  return false; // dummy return
}

/*! compare two storages, they are same if all the coefficients are travelled in the same way
    - storage pointers are the same
    - storage types are the same, sizes are the same and storage structure are the same
*/
bool sameStorage(const MatrixStorage& sto1,const MatrixStorage& sto2)
{
  if(&sto1==&sto2) return true;   //same storage pointers
  if(sto2.storageType()!=sto1.storageType() || sto2.accessType()!=sto1.accessType()) return false;
  if(sto1.sameStorage(sto2)) return true;
  return false;
}

} // end of namespace xlifepp