src/operator/Operand.hpp

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

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

/*!
  \file Operand.hpp
  \author E. Lunéville
  \since 02 mar 2012
  \date 4 may 2012

  \brief Definition of the xlifepp::Operand class

  Class xlifepp::Operand describes a couple (object, operation) where object is either an xlifepp::OperatorOnFunction,
  xlifepp::OperatorOnKernel or a xlifepp::Value (scalar, xlifepp::Vector, xlifepp::Matrix constant) and operation is an
  algebraic operator such as Fun* or *Fun or Mat*F ...

  It is used by xlifepp::OperatorOnUnknown class.
  Normally, end users has not to deal explicitly with this class

  this class stores the couple (object,operation)
    - when object is a xlifepp::Value, in a pointer to a Value
    - when object is a xlifepp::OperatorOnFunction, in a pointer to a xlifepp::OperatorOnFunction
    - when object is a xlifepp::OperatorOnKernel, in a pointer to a xlifepp::OperatorOnKernel,

  and the algebraic operation, one among :
    - standard product   (operator *)
    - inner product      (operator |)
    - cross product      (operator ^)
    - contracted product (operator %)

  Besides, this class stores possibly additionnal operations such has conjugate and transpose
  For instance, one may have the following syntaxes (F a function, V a value):
  \verbatim
  conj(F)* or *conj(F) or tran(V)^   ...
  (conj() and tran() are functions defined in Function.hpp and Value.hpp)
  \endverbatim
*/

#ifndef OPERAND_HPP
#define OPERAND_HPP

#include "config.h"
#include "space.h"
#include "utils.h"
#include "OperatorOnKernel.hpp"
#include "OperatorOnFunction.hpp"

namespace xlifepp
{

enum AlgebraicOperator {_product, _innerProduct, _crossProduct, _contractedProduct};

/*!
  \class Operand
  to deal with syntax Value AlgebraicOperation
*/
class Operand
{
protected :
  const Value* val_p;               //!< pointer to operand object (a Value)
  OperatorOnFunction* opfun_p;      //!< pointer to OperatorOnFunction object
  OperatorOnKernel* opker_p;        //!< pointer to OperatorOnKernel object
  AlgebraicOperator operation_;     //!< operand operator
  dimPair dims_;                    //!< dimension of value or function, (1,1) by default

public :
  mutable bool conjugate_;          //!< true if the operand has to be conjugated
  mutable bool transpose_;          //!< true if the operand has to be transposed

public :
  // constructors-destructor
  Operand(const OperatorOnFunction&, AlgebraicOperator); //!< construct Operand from a OperatorOnFunction
  Operand(const OperatorOnKernel&, AlgebraicOperator);   //!< construct Operand from a OperatorOnKernel
  Operand(OperatorOnKernel&, AlgebraicOperator);         //!< construct Operand from a OperatorOnKernel
  Operand(const Function&, AlgebraicOperator);           //!< construct Operand from a Function (id on Function)
  Operand(const Kernel&, AlgebraicOperator);             //!< construct Operand from a Kernel (id on Function)
  Operand(const Value&, AlgebraicOperator);              //!< construct Operand from a Value (Value is copied)
  Operand(const Value*, AlgebraicOperator);              //!< construct Operand from a Value* (Value is not copied)
  template <typename T>
  Operand(const T&, AlgebraicOperator);                  //!< construct Operand from a scalar, a Vector or a Matrix
  Operand(const Operand&);                               //!< copy construct
  Operand& operator = (const Operand&);                  //!< assign operator
  ~Operand()                                             //! destructor
  {clear();}
  void copy(const Operand&);                             //!< copy tool
  void clear();                                          //!< delete allocated pointers and reset pointers to 0

  // accessors
  AlgebraicOperator operation() const //! returns algebraic operation
  {return operation_;}
  dimPair dims() const //! returns dimensions
  {return dims_;}

  bool isFunction() const             //! true if operand with function
  {return opfun_p!=0;}
  bool isKernel() const               //! true if operand with kernel
  {return opker_p!=0;}
  bool isValue() const                //! true if operand with value
  {return val_p!=0;}
  ValueType valueType() const;        //!< return value type of value or function
  StrucType strucType() const;        //!< return structure type of value or function
  bool normalRequired() const;        //!< true if normal involved
  bool xnormalRequired() const;       //!< true if xnormal involved
  bool ynormalRequired() const;       //!< true if ynormal involved
  bool elementRequired() const;       //!< true if element is required by any function in

  void setX(const Point& P) const     //! set X point if Function involved (usefull for Kernel)
  {if(opfun_p!=0) opfun_p->setX(P);if(opker_p!=0) opker_p->setX(P);}
  void setY(const Point& P) const     //! set Y point if Function involved (usefull for Kernel)
  {if(opfun_p!=0) opfun_p->setY(P);if(opker_p!=0) opker_p->setY(P);}
//  void setN(const Vector<real_t>* np) const     //! set normal pointer if function involved (not for Kernel)
//  {if(opfun_p!=0) opfun_p->setN(np);}

  const Value& value() const;                   //!< return the operand value as object
  // const Value& value();                         //!< return the operand value as object
  Value& value();                               //!< return the operand value as object
  const OperatorOnFunction& opfunction() const; //!< return operand object as operator on function
  const OperatorOnKernel& opkernel() const;     //!< return operand object as operator on kernel
  const OperatorOnFunction* opfunctionp() const //! return operand object as pointer to operator on function
  {return opfun_p;}
  const OperatorOnKernel* opkernelp() const     //! return operand object as pointer to operator on kernel
  {return opker_p;}
  const Function& function() const;             //!< return operand object as function
  const Kernel& kernel() const;                 //!< return operand object as kernel
  const Function* functionp() const             //! return pointer to function (0 if not)
  {if(opfun_p!=0) return opfun_p->funp(); else return 0;}
  const Kernel* kernelp() const                 //! return pointer to kernel (0 if not)
  {if(opker_p!=0) return opker_p->kernelp(); else return 0;}
  void changeKernel(Kernel* ker)                //! change Kernel in operator
  {if(opker_p!=0) opker_p->changeKernel(ker);}
  bool hasExtension() const                     //! true if opfun or opker handle an extension
  {if ((opfun_p!=0 && opfun_p->ext_p!=0) || (opker_p!=0 && opker_p->ext_p!=0)) return true;
   return false;}
  const Extension* extension() const //! return extension
  {
    if(opfun_p!=0) return opfun_p->ext_p;
    if(opker_p!=0) return opker_p->ext_p;
    return 0;
  }

  template<class T> T& valueT() const; //!< return operand object as value of T type
  template<class T>
  T& value(T&, const Point&, const Vector<real_t>* np=0, const ExtensionData* =0) const;    //!< compute value function
  template<class T>
  T& value(T&, const Point&, const Point&, const Vector<real_t>* nxp=0, const Vector<real_t>* nyp=0) const;     //!< compute value kernel
  template<class T>
  T& value(T&) const;                  //!< return const value

  template<typename T, typename R>
  Vector<T> leftEval(const Vector<R>&, dimen_t&, dimen_t&, number_t) const; //!< evaluate value operand at left
  template<typename T, typename R>
  Vector<T> rightEval(const Vector<R>&, dimen_t&, dimen_t&, number_t) const; //!< evaluate value operand at right
  template<typename T, typename R>
  Vector<T> leftEval(const Point&, const Vector<R>&, dimen_t&, dimen_t&, number_t,
                     const Vector<real_t>* =0, const ExtensionData* =0) const; //!< evaluate function operand at left
  template<typename T, typename R>
  Vector<T> rightEval(const Point&, const Vector<R>&, dimen_t&, dimen_t&, number_t,
                      const Vector<real_t>* =0, const ExtensionData* =0) const; //!< evaluate function operand at right
  template<typename T, typename R>
  Vector<T> leftEval(const Point&, const Point&, const Vector<R>&, dimen_t&, dimen_t&, number_t,
                     const Vector<real_t>* =0, const Vector<real_t>* =0) const;         //!< evaluate kernel operand at left
  template<typename T, typename R>
  Vector<T> rightEval(const Point&,const Point&, const Vector<R>&, dimen_t&, dimen_t&, number_t,
                      const Vector<real_t>* =0, const Vector<real_t>* =0) const;        //!< evaluate kernel operand at right

  void print(std::ostream&) const; //!< print Operand attributes
  void print(PrintStream& os) const {print(os.currentStream());} //!< print Operand attributes
  void printsymbolic(std::ostream&) const; //!< print Operand in symbolic form
  void printsymbolic(PrintStream& os) const {printsymbolic(os.currentStream());} //!< print Operand in symbolic form
  string_t asString() const; //!< return as symbolic string
};

// external Operand functions
bool operator==(const Operand&,const Operand&); //!< compare Operands
bool operator!=(const Operand&,const Operand&); //!< compare Operands
std::ostream& operator<<(std::ostream&, const Operand&); //!< outputs Operand attributes

// implementation of template member functions
// constructor from a scalar, a vector or a matrix, ...
template <typename T>
Operand::Operand(const T& v, AlgebraicOperator aop)
{
  operation_ = aop;
  conjugate_ = false;
  transpose_ = false;
  opfun_p = 0;
  opker_p = 0;
  // check type of v
  Value::checkTypeInList(v);
  val_p = new Value(v);
}

// return the value in a operand
template<typename T> inline
T& Operand::valueT() const
{
  if(val_p == 0)
  {
    where("Operand::valueT()");
    error("operand_notavalue");
  }
  return val_p->value<T>();
}

// compute value function
// T has to be compliant with function return type, no check !!!
template<typename T> inline
T& Operand::value(T& val, const Point& pt, const Vector<real_t>* np, const ExtensionData* extdata) const
{
  if(opfun_p!=0) return opfun_p->eval(pt, val, np, extdata);
  if(opker_p!=0)
  {
      if(opker_p->kernelp()->xpar) return opker_p->eval(pt, val, 0, np);
      else return opker_p->eval(pt, val, np, 0);
  }
  return val; // dummy return
}

// compute value kernel
// T has to be compliant with kernel return type, no check !!!
template<typename T> inline
T& Operand::value(T& val, const Point& p, const Point& q, const Vector<real_t>* nxp, const Vector<real_t>* nyp) const
{
  return opker_p->eval(p, q, val, nxp, nyp);
}

// return const value,
// T has to be compliant with value,no check !!!
template<typename T> inline
T& Operand::value(T& val) const
{
  val = val_p->value<T>();
  return val;
}

//@{
//!for template compilation reasons, fake definition of dot and crossproduct in 1D
real_t dot(const real_t& x, const real_t& y);
complex_t dot(const complex_t& x, const complex_t& y);
complex_t dot(const complex_t& x, const real_t& y);
complex_t dot(const real_t& x, const complex_t& y);
real_t crossProduct(const real_t& x, const real_t& y);
complex_t crossProduct(const complex_t& x, const complex_t& y);
complex_t crossProduct(const complex_t& x, const real_t& y);
complex_t crossProduct(const real_t& x, const complex_t& y);
//@}

/*-------------------------------------------------------------------------------------
                                   evaluate Operand
  -------------------------------------------------------------------------------------*/

//@{
//! Product by a scalar
template<typename T, typename R> inline
void evalScalarProduct(const T& val, const Vector<R>& v, Vector<T>& res)
{
  // val is a scalar and v a sequence of anything
  typename std::vector<R>::const_iterator itv=v.begin();
  res.resize(v.size());
  typename std::vector<T>::iterator itr=res.begin();
  for(; itv != v.end(); itv++, itr++) *itr = val * * itv;
}

template<>
inline void evalScalarProduct(const real_t& val, const Vector<complex_t>& v, Vector<real_t>& res)
{
  error("not_handled","evalScalarProduct(Real, Vector<Complex>, Vector<Real>)");
}

template<typename T, typename R> inline
void evalScalarProduct(const Vector<T>& vec, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  // val is a vector or a matrix and v a sequence of scalar
  d=vec.size();
  res.resize(m*d);
  n=1;
  typename std::vector<R>::const_iterator itv=v.begin();
  typename std::vector<T>::iterator itr=res.begin();
  typename std::vector<T>::const_iterator itvec;
  for(number_t k = 0; k < m; k++, itv ++)
  {
    itvec=vec.begin();
    for(number_t i=0; i<d; i++, itr++, itvec++) *itr = *itv * *itvec;
  }
}

template<>
inline void evalScalarProduct(const Vector<real_t>& vec, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalScalarProduct(Vector<Real>, Vector<Complex>, Dimen, Dimen, Number, Vector<Real>)");
}

template<typename T, typename R> inline
void evalScalarProduct(const Matrix<T>& mat, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  // val is a vector or a matrix and v a sequence of scalar
  d=mat.size();
  res.resize(m*d);
  n=mat.numberOfColumns();
  typename std::vector<R>::const_iterator itv=v.begin();
  typename std::vector<T>::iterator itr=res.begin();
  typename std::vector<T>::const_iterator itmat;
  for(number_t k = 0; k < m; k++, itv ++)
  {
    itmat=mat.begin();
    for(number_t i=0; i<d; i++, itr++, itmat++) *itr = *itv * *itmat;
  }
}

template<>
inline void evalScalarProduct(const Matrix<real_t>& mat, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalScalarProduct(Matrix<Real>, Vector<Complex>, Dimen, Dimen, Number, Vector<Real>)");
}
//@}

//@{
//! Inner product
template<typename T, typename R>
void evalInnerProduct(const Vector<T>& vec, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  // vec is a vector and v a sequence of vectors of size > val.size()
  typename std::vector<R>::const_iterator itv=v.begin();
  res.resize(m);
  typename std::vector<T>::iterator itr=res.begin();
  for(number_t k = 0; k < m; k++, itr++, itv += d)
  {
    T zero(0);
    *itr = std::inner_product(vec.begin(), vec.end(), itv, zero);
  }
  d=1;
  n=1;  //new block size
}

template<>
inline void evalInnerProduct(const Vector<real_t>& vec, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalInnerProduct(Vector<Real>, Vector<Complex>, Dimen, Dimen, Number, Vector<Real>)");
}
//@}

//@{
//! Cross product
template<typename T, typename R>
void evalCrossProduct(const Vector<T>& vec, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res, bool right)
{
  // vec is a vector and v a sequence of vectors of size > val.size()
  typename std::vector<R>::const_iterator itv;
  typename std::vector<T>::iterator itr, itr2;
  dimen_t s = 1;
  if(d == 3) s = 3; //3D
  res.resize(m * s);
  itr = res.begin();
  itv = v.begin();
  for(number_t k = 0; k < m; k++, itv += d)
  {
    itr2=itr;
    crossProduct(vec, itv, itr);
    itr++;   //the last is not incremented in crossProduct
    if(right)
      for(number_t j = 0; j < s; j++, itr2++) *itr2 *= -1; //reverse sign v x u = -u x v
  }
  d=s;
  n=1; //new block size
}

template<>
inline void evalCrossProduct(const Vector<real_t>& vec, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res, bool right)
{
  error("not_handled","evalCrossProduct(Vector<Real>, Vector<Complex>, Dimen, Dimen, Number, Vector<Real>)");
}
//@}

//@{
//! Matrix product
template<typename T, typename R> inline
void evalMatrixVectorProduct(const Matrix<T>& mat, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  // mat is a matrix nbr x nbc and v a sequence of m vectors of size d >=  nbc (n should be 1)
  dimen_t s = mat.numberOfRows();
  res.resize(m * s);
  typename std::vector<R>::const_iterator itv=v.begin();
  typename std::vector<T>::iterator itr=res.begin();
  for(number_t k = 0; k < m; k++, itr += s, itv += d)
      matvec(mat.begin(), itv, itv + d, itr, itr + s);
  d=s; n=1; //new block size
}

template<>
inline void evalMatrixVectorProduct(const Matrix<real_t>& mat, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalMatrixVectorProduct(Matrix<Real>, Vector<Complex, Dimen, Dimen, Number, Vector<Real>)");
}

template<typename T, typename R> inline
void evalVectorMatrixProduct(const Matrix<T>& mat, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  // mat is a matrix and v a sequence of vectors of size >= number of matrix rows
  dimen_t s = mat.numberOfColumns();
  res.resize(m * s);
  typename std::vector<R>::const_iterator itv=v.begin();
  typename std::vector<T>::iterator itr=res.begin();
  for(number_t k = 0; k < m; k++, itr += s, itv += d)
      vecmat(mat.begin(), itv, itv + d, itr, itr + s);
  d=s; n=1; //new block size
}

template<>
inline void evalVectorMatrixProduct(const Matrix<real_t>& mat, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalMatrixProduct(Matrix<Real>, Vector<Complex>, Dimen, Dimen, Number, Vector<Real>)");
}
//@}

//@{
//! Matrix product using vector representation of matrix
template<typename T, typename R> inline
void evalMatrixVectorProduct(const Vector<T>& vec, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  //assume v is a matrices sequence (d=n*n) and vec a vector of size >= n
  res.resize(m * n);
  typename std::vector<R>::const_iterator itv=v.begin();
  typename std::vector<T>::const_iterator itvec=vec.begin();
  typename std::vector<T>::iterator itr=res.begin();
  for(number_t k = 0; k < m; k++, itr += n, itv += d)
    matvec(itv, itvec, itvec + n, itr, itr + n);
  d=n;
  n=1; //new block size
}

template<>
inline void evalMatrixVectorProduct(const Vector<real_t>& vec, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalMatrixVectorProduct(Vector<Real>, Vector<COmplex>, Dimen, Dimen, Number, Vector<Real>)");
}

template<typename T, typename R> inline
void evalVectorMatrixProduct(const Vector<T>& vec, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  //assume v is a matrices sequence (d=n*n) and vec a vector of size >= n
  res.resize(m * n);
  typename std::vector<R>::const_iterator itv=v.begin();
  typename std::vector<T>::const_iterator itvec=vec.begin();
  typename std::vector<T>::iterator itr=res.begin();
  for(number_t k = 0; k < m; k++, itr += n, itv += d)
    vecmat(itv, itvec, itvec + n, itr, itr + n);
  d=n;
  n=1; //new block size
}

template<>
inline void evalVectorMatrixProduct(const Vector<real_t>& vec, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalVectorMatrixProduct(Vector<Real>, Vector<Complex>, Dimen, Dimen, Number, Vector<Real>)");
}

template<typename T, typename R> inline
void evalMatrixMatrixProduct(const Matrix<T>& mat, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  // product mat * v_m : mat is a matrix and v_m a matrice d/n x n
  dimen_t nbr = mat.numberOfRows(), nbc=mat.numberOfColumns(), s=nbr*n;
  res.resize(m * s);
  typename std::vector<R>::const_iterator itv=v.begin();
  typename std::vector<T>::iterator itr=res.begin();
  typename Matrix<T>::const_iterator itm = mat.begin();
  for(number_t k = 0; k < m; k++, itr += s, itv += d)
    matmat(itm, nbc, itv, nbr, n, itr);
  d=s; //new block size
}

template<>
inline void evalMatrixMatrixProduct(const Matrix<real_t>& vec, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalMatrixMatrixProduct(Matrix<Real>, Vector<Complex, Dimen, Dimen, Number, Vector<Real>)");
}

template<typename T, typename R> inline
void evalMatrixMatrixProduct2(const Matrix<T>& mat, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  // product v_m * mat: mat is a matrix and v_m a matrice d/n x n
  dimen_t nbr = d/n, nbc=mat.numberOfColumns(), s=nbr*nbc;
  res.resize(m * s);
  typename std::vector<R>::const_iterator itv=v.begin();
  typename std::vector<T>::iterator itr=res.begin();
  typename Matrix<T>::const_iterator itm = mat.begin();
  for(number_t k = 0; k < m; k++, itr += s, itv += d)
    matmat(itv, n, itm, nbr, nbc, itr);
  d=s;
  n=nbr; //new block size
}

template<>
inline void evalMatrixMatrixProduct2(const Matrix<real_t>& mat, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalMatrixMatrixProduct2(Matrix<Real, Vector<Complex>, Dimen, Dimen, Number, Vector<Real>)");
}
//@}

//@{
//! Contracted product
template<typename T, typename R> inline
void evalContractedProduct(const Matrix<T>& mat, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m, Vector<T>& res)
{
  /* it is assumed that vector entry v is a sequence of m dense SQUARE matrices of size n (row storage and n*n=d)
     we consider 3 cases (for each matrix)
        - operand is a matrix nxn :  standard contracted product returning a scalar r = sum_ij A_ij*E_ij = trace(A*B)
        - operand is a n(n+1)/2 square matrix : symmetric extended contracting product returning a symmetric nxn matrix, say S
          n=3 :  [S11 S22 S33 S12 S13 S23] = A *[E11 E22 E33 E32 E31 E21]        (Voigt notation)
          n=2 :  [S11 S22 S12] = A *[E11 E22 E12]
          n=1 :  [S11] = A *[E11]
        - operand is a matrix of n^2 x n^2 : full extended contracting product returning a nxn matrix
          Sij =  C[(i-1)n+1:in;(j-1)n+1:jn]%E (block contracted product)
     the last 2 cases are generally involved in anisotropic elasticity computation, case 2 corresponds to the Voigt notation
  */

  typename std::vector<R>::const_iterator itv;
  typename std::vector<T>::iterator itr;
  number_t nA=mat.numberOfRows();

  if(nA==1)  // same as scalar * E
  {
    T val= mat[0];
    res.resize(m);
    itr=res.begin();
    for(itv = v.begin(); itv != v.end(); itv++, itr++)  *itr = val * *itv;
    return;
  }
  if(nA == n) //standard contracted product
  {
    res.resize(m);
    itv = v.begin();
    itr=res.begin();
    for(number_t k = 0; k < m; k++, itr++, itv += d)
      {
        T zero(0);
        *itr = std::inner_product(mat.begin(), mat.end(), itv, zero);
      }
    d=1;
    n=1; //new block size
    return;
  }

  if(nA == number_t(n*n) ) //full extended contracted product
  {
    res.resize(m*d);
    itv = v.begin();
    itr=res.begin();
    number_t n2=n*n;
    typename std::vector<R>::const_iterator it;
    typename Vector<T>::const_iterator itAb= mat.begin(), itA, itAij;
    for(number_t k = 0; k < m; k++, itv += d)  //block m
      for(number_t i=0; i<n; i++)
        for(number_t j=0; j<n; j++, itr++)   //compute Sij
        {
          *itr=T(0);
          it=itv;
          itAij=itAb + (i*d+j*n);
          for(number_t p=0; p<n; p++)      //contracted product
          {
            itA =itAij +p*n2;
            for(number_t q=0; q<n; q++, it++, itA++)
              *itr += *itA** it;
          }
        }
    //unchanged block size
    return;
  }

  if(2*nA == number_t(n*(n+1))) //symmetric extended contracted product
  {
    res.resize(m*d);
    itv = v.begin();
    itr=res.begin();
    Vector<R> E(nA);
    Vector<T> S(nA);
    typename Vector<R>::iterator itE;
    typename Vector<T>::iterator itS;
    for(number_t k = 0; k < m; k++, itv += d)  //block k
    {
      if(n==2)
      {
        itE=E.begin();
        *itE = *itv; itE++;     //E11
        *itE = *(itv+3); itE++; //E22
        *itE = *(itv+1);        //E21
        matvec(mat.begin(), E.begin(), E.end(), S.begin(), S.end());
        itS=S.begin();
        *itr= *itS;     itr++; //S11
        *itr= *(itS+2); itr++; //S12
        *itr= *(itS+2); itr++; //S21
        *itr= *(itS+1); itr++; //S22
      }
      else
      {
        itE=E.begin();
        *itE = *itv; itE++;     //E11
        *itE = *(itv+4); itE++; //E22
        *itE = *(itv+8); itE++; //E33
        *itE = *(itv+5); itE++; //E23
        *itE = *(itv+2); itE++; //E13
        *itE = *(itv+1);        //E12
        matvec(mat.begin(), E.begin(), E.end(), S.begin(), S.end());
        itS=S.begin();
        *itr= *itS;      itr++; //S11
        *itr = *(itS+5); itr++; //S12
        *itr = *(itS+4); itr++; //S13
        *itr = *(itS+5); itr++; //S21
        *itr = *(itS+1); itr++; //S22
        *itr = *(itS+3); itr++; //S23
        *itr = *(itS+4); itr++; //S31
        *itr = *(itS+3); itr++; //S32
        *itr = *(itS+2); itr++; //S33
      }
    }
    //unchanged block size
    return;
  }
}

template<>
inline void evalContractedProduct(const Matrix<real_t>& mat, const Vector<complex_t>& v, dimen_t& d, dimen_t& n, number_t m, Vector<real_t>& res)
{
  error("not_handled","evalContractedProduct(Matrix<Real>, Vector<Complex>, Dimen, Dimen, Number, Vector<Real>)");
}
//@}

//-----------------------------------------------------
// left and right operand action
//-----------------------------------------------------
//@{
/*! evaluate operand with value or function at left or right applied to a collection of scalar/vector/matrix values of m basis functions
    (v_1)_1,...,(v_1)_d,(v_2)_1,...,(v_2)_d, ..., (v_m)_1,...,(v_m)_d
    return in same structure w= op aop v
    (w_1)_1,...,(w_1)_q,(w_2)_1,...,(w_2)_q, ..., (w_m)_1,...,(w_m)_q

    v  : basis function values
    m  : number of blocks (each block corresponds to a basis function)
    d  : size of block (may be updated)
    n  : dimension of matrix, when blocks are matrices (p=d/n):
         [v_1,...,v_d] = [v_11,...,v_1n, v_21,...,v_2n, ..., v_p1,...,v_pn]

    Note : no consistency test are done here !!!
    Note : d and n may be updated by operations
*/
template<typename T, typename R> inline
Vector<T> Operand::leftEval(const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m) const //! value op V
{
  Vector<T> r;
  StrucType sto=strucType();  //struct of operand (scalar, vector, matrix)
  ValueType vt=valueType();
  if(sto == _scalar) //scalar value in operand all operations degenerate as a scalar product
  {
    T val=val_p->value<T>();
    if(vt==_complex && conjugate_) val=conj(val);
    evalScalarProduct(val, v, r);
    return r;
  }

  if(sto == _vector)
  {
    Vector<T> vec= val_p->value<Vector<T> >();
    if(vt==_complex && conjugate_) vec=conj(vec);
    switch(operation_)
    {
      case _innerProduct : //not hermitian product
      {
        evalInnerProduct(vec, v, d, n, m, r);
        return r;
      }
      case _crossProduct :  //in 3D or 2D (return a scalar)
      {
        evalCrossProduct(vec, v, d, n, m, r, false);
        return r;
      }
      case _product : //vector * scalar or vector * matrix
      {
        if (d==1) //vector * scalar
        {
          evalScalarProduct(vec, v, d, n, m, r);
          return r;
        }
        else if (n>=1) //vector * matrix
        {
          evalVectorMatrixProduct(vec, v, d, n, m, r);
          return r;
        }
      }
      default:
        break;
    }
  }

  if(sto == _matrix)
  {
    Matrix<T> mat= val_p->value<Matrix<T> >();
    if(vt==_complex && conjugate_) mat=conj(mat);
    if(vt==_real && transpose_) mat.transpose();
    switch(operation_)
    {
      case _product :
      {
        if (d==1) //vector * scalar
        {
          evalScalarProduct(mat, v, d, n, m, r);
          return r;
        }
        else
        {
          if (n==1) { evalMatrixVectorProduct(mat, v, d, n, m, r); return r; }
          else { evalMatrixMatrixProduct(mat, v, d, n, m, r); return r; }
        }
      }
      case _contractedProduct:
      {
        evalContractedProduct(mat, v, d, n, m, r);
        return r;
      }
      default:
        break;
    }
  }

  error("not_handled", "Operand::leftEval");
  return r;
}

template<typename T, typename R> inline
Vector<T> Operand::rightEval(const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m) const //! V op value
{
  Vector<T> r;
  StrucType sto=strucType();  //struct of operand (scalar, vector, matrix)
  ValueType vt=valueType();
  if (sto == _scalar) //scalar value in operand all operations degenerate as a scalar product
  {
    T val=val_p->value<T>();
    if(vt==_complex && conjugate_) val=conj(val);
    evalScalarProduct(val, v, r);
    return r;
  }

  if (sto == _vector)
  {
    Vector<T> vec= val_p->value<Vector<T> >();
    if(vt==_complex && conjugate_) vec=conj(vec);
    switch (operation_)
    {
      case _innerProduct : //not hermitian product
      {
        evalInnerProduct(vec, v, d, n, m, r);
        return r;
      }
      case _crossProduct :  //in 3D or 2D (return a scalar)
      {
        evalCrossProduct(vec, v, d, n, m, r, true);
        return r;
      }
      case _product : //vector * scalar or vector * matrix
      {
        if (d==1) //vector * scalar
        {
          evalScalarProduct(vec, v, d, n, m, r);
          return r;
        }
        else if (n>=1) //vector * matrix
        {
          evalMatrixVectorProduct(vec, v, d, n, m, r);
          return r;
        }
      }
      default:
        break;
    }
  }

  if (sto == _matrix)
  {
    Matrix<T> mat= val_p->value<Matrix<T> >();
    if (vt==_complex && conjugate_) mat=conj(mat);
    if (vt==_real && transpose_) mat.transpose();
    switch (operation_)
    {
      case _product :
      {
        if (d==1) //vector * scalar
        {
          evalScalarProduct(mat, v, d, n, m, r);
          return r;
        }
        else
        {
          if (n==1) { evalVectorMatrixProduct(mat, v, d, n, m, r); return r; }
          else { evalMatrixMatrixProduct2(mat, v, d, n, m, r); return r; }
        }
      }
      case _contractedProduct :
      {
        evalContractedProduct(mat, v, d, n, m, r);
        return r;
      }
      default:
        break;
    }
  }

  error("not_handled", "Operand::leftEval");
  return r;
}

//=========================================================================================
//                               operand with function
//=========================================================================================

template<typename T, typename R> inline
Vector<T> Operand::leftEval(const Point& pt, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m,
                            const Vector<real_t>* np, const ExtensionData* extdata) const //! f(p) op V
{
  Vector<T> r;
  StrucType sto=strucType();  //struct of operand (scalar, vector, matrix)
  ValueType vt=valueType();
  if (sto == _scalar) //scalar value in operand or vector, all operations degenerate as a scalar product
  {
    T val;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       real_t valR;
       value(valR, pt, np, extdata);
       val=valR;
    }
    else value(val, pt, np, extdata); //result is mandatory complex
    if (vt==_complex && conjugate_) val=conj(val);
    evalScalarProduct(val,v,r);
    return r;
  }

  if (sto == _vector)
  {
    Vector<T> vec;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       Vector<real_t> vecR;
       value(vecR, pt, np, extdata);
       vec=vecR;
    }
    else value(vec, pt, np); //result is mandatory complex
    if (vt==_complex && conjugate_) vec=conj(vec);
    switch (operation_)
    {
      case _innerProduct : //not hermitian product
      {
        evalInnerProduct(vec, v, d, n, m, r);
        return r;
      }
      case _crossProduct :  //in 3D or 2D (return a scalar)
      {
        evalCrossProduct(vec, v, d, n, m, r, false);
        return r;
      }
      case _product : //only vector * matrix
      {
        if (d==1) //vector * scalar
        {
          evalScalarProduct(vec, v, d, n, m, r);
          return r;
        }
        else if (n>=1) //vector * matrix
        {
          evalVectorMatrixProduct(vec, v, d, n, m, r);
          return r;
        }
      }
      default :
          break;
    }
  }

  if(sto == _matrix)
  {
    Matrix<T> mat;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       Matrix<real_t> matR;
       value(matR, pt, np, extdata);
       mat=matR;
    }
    else value(mat, pt, np, extdata);//result is mandatory complex
    if (vt==_complex && conjugate_) mat=conj(mat);
    if (vt==_real && transpose_) mat.transpose();
    switch (operation_)
    {
      case _product : // matrix * vector
      {
        if (d==1) //vector * scalar
        {
          evalScalarProduct(mat, v, d, n, m, r);
          return r;
        }
        else
        {
          //if (n==mat.numberOfColumns()) { evalMatrixVectorProduct(mat, v, d, n, m, r); return r; }
          if (n==1) { evalMatrixVectorProduct(mat, v, d, n, m, r); return r; }
          else { evalMatrixMatrixProduct(mat, v, d, n, m, r); return r; }
        }
      }
      case _contractedProduct :
      {
        evalContractedProduct(mat, v, d, n, m, r);
        return r;
      }
      default:
        break;
    }
  }

  error("not_handled", "Operand::leftEval");
  return r;
}


template<typename T, typename R> inline
Vector<T> Operand::rightEval(const Point& pt, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m,
                             const Vector<real_t>* np, const ExtensionData* extdata) const //!  V op f(p)
{
  Vector<T> r;
  StrucType sto=strucType();  //struct of operand (scalar, vector, matrix)
  ValueType vt=valueType();
  if (sto == _scalar) //scalar value in operand or vector, all operations degenerate as a scalar product
  {
    T val;
    if(vt==_real)  //to deal with real matrix and complex result (T)
    {
       real_t valR;
       value(valR, pt, np, extdata);
       val=valR;
    }
    else value(val, pt, np, extdata); //result is mandatory complex
    if (vt==_complex && conjugate_) val=conj(val);
    evalScalarProduct(val,v,r);
    return r;
  }

  if (sto == _vector)
  {
    Vector<T> vec;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       Vector<real_t> vecR;
       value(vecR, pt, np, extdata);
       vec=vecR;
    }
    else value(vec, pt, np, extdata); //result is mandatory complex
    if (vt==_complex && conjugate_) vec=conj(vec);
    switch(operation_)
    {
      case _innerProduct : //not hermitian product
        evalInnerProduct(vec, v, d, n, m, r);
        return r;
      case _crossProduct :  //in 3D or 2D (return a scalar)
        evalCrossProduct(vec, v, d, n, m, r, false);
        return r;
      case _product : //only vector * matrix
        if (d==1) //vector * scalar
        { evalScalarProduct(vec, v, d, n, m, r); return r; }
        else if (n>=1) //vector * matrix
        { evalMatrixVectorProduct(vec, v, d, n, m, r); return r; }
      default:
        break;
    }
  }

  if(sto == _matrix)
  {
    Matrix<T> mat;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       Matrix<real_t> matR;
       value(matR, pt, np, extdata);
       mat=matR;
    }
    else value(mat, pt, np, extdata); //result is mandatory complex
    if (vt==_complex && conjugate_) mat=conj(mat);
    if (vt==_real && transpose_) mat.transpose();
    switch(operation_)
    {
      case _product : // matrix * vector
        if (d==1) //vector * scalar
        { evalScalarProduct(mat, v, d, n, m, r); return r; }
        else
        {
          //if (n==mat.numberOfRows()) { evalVectorMatrixProduct(mat, v, d, n, m, r); return r; }
          if (n==1) { evalVectorMatrixProduct(mat, v, d, n, m, r); return r; }
          else { evalMatrixMatrixProduct2(mat, v, d, n, m, r); return r; }
        }
      case _contractedProduct :
        evalContractedProduct(mat, v, d, n, m, r);
        return r;
      default:
        break;
    }
  }

  error("not_handled", "Operand::rightEval");
  return r;
}

//=========================================================================================
//                               operand with kernel
//=========================================================================================
//  pointer opker cannot be a null pointer !!!

template<typename T, typename R> inline
Vector<T> Operand::leftEval(const Point& x, const Point& y, const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m,
                            const Vector<real_t>* nx, const Vector<real_t>* ny) const //! f(p) op V
{
  Vector<T> r;
  StrucType sto=strucType();  //struct of operand (scalar, vector, matrix)
  ValueType vt=valueType();
  if (sto == _scalar) //scalar value in operand or vector, all operations degenerate as a scalar product
  {
    T val;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       real_t valR;
       opker_p->eval(x, y, valR, nx, ny);
       val=valR;
    }
    else opker_p->eval(x, y, val, nx, ny);  //result is mandatory complex
    if (vt==_complex && conjugate_) val=conj(val);
    evalScalarProduct(val,v,r);
    return r;
  }

  if (sto == _vector)
  {
    Vector<T> vec;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       Vector<real_t> vecR;
       opker_p->eval(x, y, vecR, nx, ny);
       vec=vecR;
    }
    else opker_p->eval(x, y, vec, nx, ny); //result is mandatory complex
    if (vt==_complex && conjugate_) vec=conj(vec);
    switch (operation_)
    {
      case _innerProduct : //not hermitian product
      {
        evalInnerProduct(vec, v, d, n, m, r);
        return r;
      }
      case _crossProduct :  //in 3D or 2D (return a scalar)
      {
        evalCrossProduct(vec, v, d, n, m, r, false);
        return r;
      }
      case _product :
      {
        if (d==1) //vector * scalar
        {
          evalScalarProduct(vec, v, d, n, m, r);
          return r;
        }
        else if (n>=1) //vector * matrix
        {
          evalVectorMatrixProduct(vec, v, d, n, m, r);
          return r;
        }
      }
      default :
          break;
    }
  }

  if(sto == _matrix) // compute kernel matrix in a vector
  {
    dimPair dims(0,0);
    Matrix<T> mat;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       Vector<real_t> matR;
       opker_p->eval(x, y, matR, nx, ny, &dims);
       mat = Matrix<T>(dims.first, dims.second, matR.begin(), matR.end());
    }
    else
       {
        Vector<T> matT;
        opker_p->eval(x, y, matT, nx, ny, &dims); //result is mandatory complex
        mat = Matrix<T>(dims.first, dims.second, matT.begin(), matT.end());
       }

    if (vt==_complex && conjugate_) mat.conjugate();
    if (transpose_) mat.transpose();
    switch (operation_)
    {
      case _product : // matrix * vector
      {
        if (d==1) //vector * scalar
        {
          evalScalarProduct(mat, v, d, n, m, r);
          return r;
        }
        else
        {
          //if (n==mat.numberOfColumns()) { evalMatrixVectorProduct(mat, v, d, n, m, r); return r; }
          if (n==1) { evalMatrixVectorProduct(mat, v, d, n, m, r); return r; }
          else { evalMatrixMatrixProduct(mat, v, d, n, m, r); return r; }
        }
      }
      case _contractedProduct :
      {
        evalContractedProduct(mat, v, d, n, m, r);
        return r;
      }
      default:
        break;
    }
  }

  error("not_handled", "Operand::leftEval");
  return r;
}


template<typename T, typename R> inline
Vector<T> Operand::rightEval(const Point& x, const Point& y,  const Vector<R>& v, dimen_t& d, dimen_t& n, number_t m,
                             const Vector<real_t>* nx, const Vector<real_t>* ny) const //!  V op f(p)
{
  Vector<T> r;
  StrucType sto=strucType();  //struct of operand (scalar, vector, matrix)
  ValueType vt=valueType();
  if (sto == _scalar) //scalar value in operand or vector, all operations degenerate as a scalar product
  {
    T val;
    if(vt==_real)  //to deal with real matrix and complex result (T)
    {
       real_t valR;
       opker_p->eval(x, y, valR, nx, ny);
       val=valR;
    }
    else opker_p->eval(x, y, val, nx, ny); //result is mandatory complex
    if (vt==_complex && conjugate_) val=conj(val);
    evalScalarProduct(val,v,r);
    return r;
  }

  if (sto == _vector)
  {
    Vector<T> vec;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       Vector<real_t> vecR;
       opker_p->eval(x, y, vecR, nx, ny);
       vec=vecR;
    }
    else opker_p->eval(x, y, vec, nx, ny); //result is mandatory complex
    if (vt==_complex && conjugate_) vec=conj(vec);
    switch(operation_)
    {
      case _innerProduct : //not hermitian product
        evalInnerProduct(vec, v, d, n, m, r);
        return r;
      case _crossProduct :  //in 3D or 2D (return a scalar)
        evalCrossProduct(vec, v, d, n, m, r, false);
        return r;
      case _product : //only vector * matrix
        if (d==1) //vector * scalar
        { evalScalarProduct(vec, v, d, n, m, r); return r; }
        else if (n>=1) //vector * matrix
        { evalMatrixVectorProduct(vec, v, d, n, m, r); return r; }
      default:
        break;
    }
  }

  if(sto == _matrix)
  {
    Matrix<T> mat;
    if(vt==_real)  //to deal with real matrix and complex result
    {
       Matrix<real_t> matR;
       opker_p->eval(x, y, matR, nx, ny);
       mat=matR;
    }
    else opker_p->eval(x, y, mat, nx, ny);//result is mandatory complex
    if (vt==_complex && conjugate_) mat=conj(mat);
    if (vt==_real && transpose_) mat.transpose();
    switch(operation_)
    {
      case _product : // matrix * vector
        if (d==1) //vector * scalar
        { evalScalarProduct(mat, v, d, n, m, r); return r; }
        else
        {
          //if (n==mat.numberOfRows()) { evalVectorMatrixProduct(mat, v, d, n, m, r); return r; }
          if (n==1) { evalVectorMatrixProduct(mat, v, d, n, m, r); return r; }
          else { evalMatrixMatrixProduct2(mat, v, d, n, m, r); return r; }
        }
      case _contractedProduct :
        evalContractedProduct(mat, v, d, n, m, r);
        return r;
      default:
        break;
    }
  }

  error("not_handled", "Operand::rightEval");
  return r;
}
//@}

/*! evaluate "tensor" operation between 2 small vectors, and add it to a given a nu x nv matrix mat (stored by rows)

    op : algebraic operator involved (either product, inner product, cross product, contracted product)
    alpha is a coefficient applied to the added "tensor" matrix
    u (resp. v) is
     - either a list of nu scalar values (resp. nv scalar values)
     - or a list of nbu (resp. nbv) vector values (u_1,u_2,..u_nu) stored as
         (u_1)1 ...(u_1)mu (u_2)1 ... (u_2)mu .....(u_nu)1 ...(u_nu)mu
          with mu=u.size()/nu
         for scalar case (mu=mv=1) all product are equivalent
         for vector case (mu=mv>1) only inner product is currently available
         for matrix case (mu=mv>1) only contracted product is currently available

*/
template <typename T, typename KU, typename KV> inline
Matrix<T>& tensorOpAdd(const AlgebraicOperator& op, const std::vector<KU>& u, number_t nu,
                      const std::vector<KV>& v, number_t nv, Matrix<T>& mat, const T& alpha)
{

  typename Matrix<T>::iterator itm = mat.begin();
  typename std::vector<KU>::const_iterator itu;
  typename std::vector<KV>::const_iterator itv;
  number_t mu = u.size() / nu, mv = v.size() / nv;
  if(mu != mv)
  {
    where("tensorOpAdd(...)");
    error("bad_size", mv, mu);
  }

  if(mu == 1) //scalar case
  {
    for(itu = u.begin(); itu != u.end(); itu++)
      for(itv = v.begin(); itv != v.end(); itv++, itm++)
        *itm += alpha** itu** itv;
  }
  else
    switch(op)
    {
    case _innerProduct :
    {
      for(itu = u.begin(); itu != u.end(); itu += mu)
        for(itv = v.begin(); itv != v.end(); itv += mv, itm++)
          {
            T zero(0);
            *itm += alpha * std::inner_product(itu, itu + mu, itv, zero);
          }
      return mat;
      break;
    }
    case _contractedProduct:
    {
      for(itu = u.begin(); itu != u.end(); itu += mu)
        for(itv = v.begin(); itv != v.end(); itv += mv, itm++)
          for(number_t i=0; i<mu; i++)
            *itm += alpha * *(itu+i) * *(itv+i);
      return mat;
      break;
    }
    case _product :
    case _crossProduct :
      error("not_handled", "tensorOpAdd(...)");
  }
  return mat; //fake return
}

//template <typename T> inline
//Matrix<T>& tensorOpAdd(const AlgebraicOperator& op, const std::vector<T>& u, number_t nu,
//                      const std::vector<T>& v, number_t nv, Matrix<T>& mat, const T& alpha)
//{
//
//  typename Matrix<T>::iterator itm = mat.begin();
//  typename std::vector<T>::const_iterator itu, itv;
//  number_t mu = dimen_t(u.size()) / nu, mv = dimen_t(v.size()) / nv;
//  if(mu != mv) error("free_error", "in tensorOpAdd(...) : not the same size of vectors");
//
//  if(mu == 1) //scalar case
//    {
//      for(itu = u.begin(); itu != u.end(); itu++)
//        for(itv = v.begin(); itv != v.end(); itv++, itm++)
//          *itm += alpha** itu** itv;
//    }
//  else
//    switch(op)
//      {
//      case _innerProduct :
//        {
//          for(itu = u.begin(); itu != u.end(); itu += mu)
//            for(itv = v.begin(); itv != v.end(); itv += mv, itm++)
//              {
//                T zero(0);
//                *itm += alpha * std::inner_product(itu, itu + mu, itv, zero);
//              }
//          return mat;
//        }
//        break;
//      case _contractedProduct:
//        {
//          for(itu = u.begin(); itu != u.end(); itu += mu)
//            for(itv = v.begin(); itv != v.end(); itv += mv, itm++)
//              for(number_t i=0; i<mu; i++)
//                *itm += alpha * *(itu+i) * *(itv+i);
//          return mat;
//        }
//        break;
//      case _product :
//      case _crossProduct :
//        error("free_error", "in tensorOp(...) : operator on vectors not available");
//      }
//  return mat; //fake return
//}


//same function with iterator instead of Matrix<T>
//template <typename T, typename iterator> inline
//void tensorOpAdd(const AlgebraicOperator& op, const std::vector<T>& u, number_t nu,
//                 const std::vector<T>& v, number_t nv, iterator itm, const T& alpha)
//{
//
//  typename std::vector<T>::const_iterator itu, itv;
//  number_t mu = dimen_t(u.size()) / nu, mv = dimen_t(v.size()) / nv;
//  if(mu != mv) error("free_error", "in tensorOpAdd(...) : not the same size of vectors");
//  if(mu == 1) //scalar case
//    {
//      for(itu = u.begin(); itu != u.end(); itu++)
//        for(itv = v.begin(); itv != v.end(); itv++, itm++)
//          *itm += alpha** itu** itv;
//    }
//  else
//    switch(op)
//      {
//      case _innerProduct :
//        {
//          for(itu = u.begin(); itu != u.end(); itu += mu)
//            for(itv = v.begin(); itv != v.end(); itv += mv, itm++)
//              {
//                T zero(0);
//                *itm += alpha * std::inner_product(itu, itu + mu, itv, zero);
//              }
//        }
//        break;
//      case _product :
//      case _crossProduct :
//      case _contractedProduct:
//        error("free_error", "in tensorOp(...) : operator on vectors not available");
//      }
//}

/*! evaluate "tensor" operation between 2 small vectors and "coefficient" ker, and add it to a given a nu x nv matrix mat (stored by rows)
                                   u opu ker opv  v
    u (resp. v) is
     - either a list of nu scalar values (resp. nv scalar values)
     - or a list of nbu (resp. nbv) vector values (u_1,u_2,..u_nu) stored as
         (u_1)1 ...(u_1)mu (u_2)1 ... (u_2)mu .....(u_nu)1 ...(u_nu)mu
          with mu=u.size()/nu
    for scalar case (mu=mv=1) all product are equivalent
    for vector case (mu=mv>1) only inner product is currently available

    itk  is  an iterator to the kernel part (either a scalar, a vector (mk or nk) or a matrix (mk x nk))
*/
template <typename T, typename iteratorM, typename iteratorK> inline
void tensorOpAdd(const AlgebraicOperator& opu, const AlgebraicOperator& opv, const std::vector<T>& u, number_t nu,
                 const std::vector<T>& v, number_t nv, iteratorM itm, iteratorK itk, number_t mk, number_t nk)
{

  typename std::vector<T>::const_iterator itu, itv;
  number_t mu = dimen_t(u.size()) / nu, mv = dimen_t(v.size()) / nv;
  if(mu != mv)
  {
    where("tensorOpAdd(...)");
    error("bad_size", mv, mu);
  }
  if(mu == 1) //scalar case, all, operations are product
  {
    for(itu = u.begin(); itu != u.end(); ++itu)
      for(itv = v.begin(); itv != v.end(); ++itv, ++itm)
        *itm +=  *itu * *itk * *itv;
    return;
  }

  /*vector case mu = mv > 1,  supported operations :
        U | alpha * V  means   U | (alpha * V)   alpha a scalar or a matrix mu x mv
        U * alpha | V  means  (U * alpha) | V    alpha a scalar or a matrix mu x mv
        U | alpha ^ V  means   U | (alpha ^ V)   alpha a vector mv x 1
        U ^ alpha | V  means  (U ^ alpha) | V    alpha a vector mu x 1
  */
  if(opu==_innerProduct && opv==_product)
  {
    for(itu = u.begin(); itu != u.end(); itu += mu)
      for(itv = v.begin(); itv != v.end(); itv += mv, itm++)
      {
        // T zero(0);
        //*itm += alpha * std::inner_product(itu, itu + mu, itv, zero);
      }
    return;
  }

  if(opu==_innerProduct && opv==_product)
  {
    return;
  }

  if(opu==_innerProduct && opv==_product)
  {
    return;
  }

  error("not_handled", "tensorOpAdd(...)");
}

} // end of namespace xlifepp

#endif /* OPERAND_HPP */