contrib/bamg/RNM.hpp

// ********** DO NOT REMOVE THIS BANNER **********
// ORIG-DATE:    29 fev 2000
// -*- Mode : c++ -*-
//
// SUMMARY  : array modelisation
// USAGE    : LGPL
// ORG      : LJLL Universite Pierre et Marie Curie, Paris,  FRANCE
// AUTHOR   : Frederic Hecht
// E-MAIL   : frederic.hecht@ann.jussieu.fr
//

/*


 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.

 Freefem++  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 Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA


 */
#ifndef KNM_H_
#define KNM_H_
// version march  2010 FH.
// add iterator for compatibility with gmm++
// ----------------------
// une tentative qui ne marche pas
// de tableau constant
#include <complex>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <cassert>


using namespace std;
#define const_R   R

#include <cstdlib>
inline void Check_Kn(const char * str,const char * file,int line)
{
 cerr << "CHECK_KN: " << str << " in file: " << file << ", line " << line <<endl;
 assert(0);
 abort();
}

#define K_bigassert(i)  if (!(i)) Check_Kn(#i,__FILE__,__LINE__);
#define RNM_FATAL_ERROR(i) Check_Kn(i,__FILE__,__LINE__);
#ifdef CHECK_KN

#define K_throwassert(i)  if (!(i)) Check_Kn(#i,__FILE__,__LINE__);

#else
#define K_throwassert(i)
#endif
// version du 29 fev 2000
//  correction   for (... lj++,ui++  qui apelle le produit scalaire
//  petite correction  throwassert
// ajoute de operateur /= et *= sur des vecteurs
//   suppression de constructeur qui pose de probleme
//   correction oper +=  ...  dans RNM_op.h ligne 56  change = en oper
// version de 25 nov 99  sans const R
// ajoute de '.' pour extraire une colonne, ou ligne  , ...
//  version du 22 nov 1999   cast to KN_<const R>
//   version du 21 nov 1999  correction delete
//  version du 13 nov 1999
//  version du 18 mars 99
//  F. Hecht
// attention les indexations les indexations peuvent changer
//  puisque que l'on peut prendre la transposer d'une matrice
// tableau
// mais ils partent de 0
// version corrigee du 15/11/98
// version avec sous tableau  ---  mars 99
// -------
//  remarque du 8 mars 99 FH
// class pour prendre des sous-tableau
// attention aux PB de continute dans les tableaux
// on a supposer que les tableaux multi indices pouvait est vue comme
// un tableau continue ce qui est generalement faux quand l'on en
// prend un sous tableau
//   exemple: un tableau 3,5 est numerote comme:
//    0  3  6  9 12
//    1  4  7 10 13
//    2  5  8 11 14
//             step
//   indexi  n 1
//   indexj  m n
//   est le sous tableau  3,3  n'est pas numeroter consecutivement
//
//    Donc la fonction  IsVector1() nous dit si un tableau
//    a un 2 ou 3 indices est ou non consecutif en memoire
//
//  --  ajoute d'une classe VirtualMatrice
// pour modeliser le produit matrice vecteur
//  x = A*v; via des fonctions virtuelle
//  ----------------------------------
//   version du 6 mars 2001 FH
//   ---  initialisation --
//   --------------------------------
//   version du 9 oct 2001 FH
//   ajoute de constructeur par defaut d'une vecteur
//   +  set , pour definir le vecteur
//   ou l'affectation (bof bof)
// ---------------------
//  version sep 2002
//  ajoute  operateur >> pour  KN<R> et KN_<R>
//  --------------------
//  version  april 2003
//  ajoute un gestion auto de
//  la fonction InternalError pour les matriceVirtuel
//  --------------------
//   version jan 2004
//   correction pour go ++
//   des operateur  #=  pour les matrices et tenseurs
//  ----------------------
//   version  feb 2004
//   v(i(:)) =  w   //  i(1:10)
//   w=u(i(:))  //
//   version mars 2004 make small correction
//   in  ITAB operator problem if non type R a defi
//   -------------------
//   Modif pour version avec les Complex   mai 2004
//   (u,v)  donne le produit complex utiliser dans le produit matrice vecteur
//   (u,conj(v))  donne le produit hermitiene pour le gradient conjugue
//
//   -- de fonction dans le cas real
// modif for g++ 4.0 and xlc++   mai 2005
//  adding some this->
//   mars 2007
// correction in operator operation:b -1*c
// aout 2007,
//  correct y = A*x ; when y is unset
//  correct y += A*x ; when y is unset
//  re-correct += sep 2007
//  add size of the matrix in VirtualMatrix class.
//   mars 2010 add  unset KNM case ...
// ----------------
inline double  conj(const double & x){return x;}
inline float  conj(const float &x){return x;}
inline long  conj(const long &x){return x;}
inline double  real(const double &x){return x;}
inline float  real(const float &x){return x;}

namespace RNM
{
  template<class T> inline T Min (const T &a,const T &b){return a < b ? a : b;}
  template<class T> inline T Max (const T &a,const T & b){return a > b ? a : b;}
  template<class T> inline T Abs (const T &a){return a <0 ? -a : a;}

  template<class T> inline void Exchange (T& a,T& b) {T c=a;a=b;b=c;}
  template<class T> inline T Max (const T &a,const T & b,const T & c){return Max(Max(a,b),c);}
  template<class T> inline T Min (const T &a,const T & b,const T & c){return Min(Min(a,b),c);}
  // specialisation cas complex ---
  template<class T>
  inline complex<T> Min(const complex<T> &a,complex<T> &b)
  { return complex<T>(min(a.real(),b.real()),min(a.imag(),b.imag()));}
  template<class T>
  inline complex<T> Max(const complex<T> &a,const complex<T> &b)
  { return complex<T>(max(a.real(),b.real()),max(a.imag(),b.imag()));}

  /*inline complex<double> Min(const complex<double> &a,complex<double> &b)
    { return complex<double>(Min(real(a),real(b)),Min(imag(a),imag(b)));}
    inline complex<double> Max(const complex<double> &a,const complex<double> &b)
    { return complex<double>(Max(real(a),real(b)),Max(imag(a),imag(b)));}
  */
}
//  ----

template<class R> class KNMK_ ;
template<class R> class KNM_ ;
template<class R> class KN_ ;
template<class R> class KN__iterator ;
template<class R> class KN__const_iterator ;
template<class R> class TKN_ ; // KN_ transpose
template<class R> class notKN_ ; // KN_ not
template<class R> class notnotKN_ ; // KN_ not not

template<class R> class KNMK ;
template<class R> class KNM ;
template<class R> class KN ;

template<class R> class conj_KN_ ;
template<class R> class Add_KN_;
template<class R> class Sub_KN_;
template<class R> class Mulc_KN_;
template<class R> class Add_Mulc_KN_;
template<class R> class Mul_KNM_KN_;
template<class R> class DotStar_KN_;
template<class R> class DotSlash_KN_;

template<class R> class outProduct_KN_;
template<class R> class if_KN_;
template<class R> class if_arth_KN_;
template<class R> class ifnot_KN_;
template<class R,class I> class KN_ITAB;

template<class R,typename A,typename B> class F_KN_;


#ifndef ffassert
#define ffassert assert
#endif

// gestion des erreur interne --
#ifndef InternalError
typedef void (* TypeofInternalErrorRoutine)(const char *) ;
static TypeofInternalErrorRoutine &InternalErrorRoutinePtr()
{
  static TypeofInternalErrorRoutine routine=0;
  return routine;
}

static void InternalError(const char * str) {
  if (InternalErrorRoutinePtr() ) (*InternalErrorRoutinePtr())(str);
  cerr << str;
  exit(1);
}
inline void  SetInternalErrorRoutine(TypeofInternalErrorRoutine f)
{
  InternalErrorRoutinePtr()=f;
}
#endif
//  --
template<class P,class Q>
   struct PplusQ { const P & p;const Q & q;
    PplusQ(const P & pp,const Q & qq) : p(pp),q(qq){}
    };

template<class R>
struct  VirtualMatrice { public:
    int N,M;
    VirtualMatrice(int nn,int mm): N(nn),M(mm) {}
    VirtualMatrice(int nn): N(nn),M(nn) {}
  //  y += A x
  virtual void addMatMul(const KN_<R> &  x, KN_<R> & y) const =0;
  virtual void addMatTransMul(const KN_<R> &  , KN_<R> & ) const
    { InternalError("VirtualMatrice::addMatTransMul not implemented "); }
  virtual bool WithSolver() const {return false;} // by default no solver
  virtual void Solve( KN_<R> &  ,const KN_<R> & ) const
    { InternalError("VirtualMatrice::solve not implemented "); }

#ifdef VersionFreeFempp
  virtual bool ChecknbLine  (int n) const= 0;
  virtual bool ChecknbColumn  (int m) const =0;
#else
  virtual bool ChecknbLine  (int n) const {return true;}
  virtual bool ChecknbColumn  (int m) const {return true;}
#endif
  struct  plusAx { const VirtualMatrice * A; const KN_<R>   x;
   plusAx( const VirtualMatrice * B,const KN_<R> &  y) :A(B),x(y)
      { ffassert(B->ChecknbColumn(y.N())); }
    };

   plusAx operator*(const KN_<R> &  x) const {return plusAx(this,x);}

  struct  plusAtx { const VirtualMatrice * A; const KN_<R>   x;
   plusAtx( const VirtualMatrice * B,const KN_<R> &  y) :A(B),x(y)
    {ffassert(B->ChecknbLine(y.N()));} };

  struct  solveAxeqb { const VirtualMatrice * A; const KN_<R>   b;
   solveAxeqb( const VirtualMatrice * B,const KN_<R> &  y) :A(B),b(y)
    {ffassert(B->ChecknbColumn(y.N()));} };

  virtual ~VirtualMatrice(){}
};


//template <class R> class MatriceCreuseMulKN_;
//template <class R> class MatriceCreuseDivKN_;

class ShapeOfArray;

class FromTo{ public:
  long from,to;
  FromTo(long i,long j):from(i),to(j) {K_throwassert(i<j);}
 };

class SubArray{ public:
  const long n,step,start;
//  SubArray(char  nn): n(-1),step(1),start(0) {}
 explicit SubArray(long nn,long sta=0,long s=1): n(nn),step(s),start(sta) {}
  SubArray(const FromTo& ft) : n(ft.to-ft.from+1),step(1),start(ft.from) {}
  SubArray(const ShapeOfArray & ); // all
  long end()  const  { return start+ step*n;}
  long last() const  { return start+ step*(n-1);}
  long len1() const  { return step*(n-1);}
};


class ShapeOfArray{ protected:
  public:

    long n;     //   n  nb of item
    long step;  //   step  nb of between 2 item
    long next;  //  the   next array of same type in matrix for subarray
              // by default  no next ( in case of KN, and KNM  -next is
              // a counter of destruction  (use in frefem++)
  ShapeOfArray(const ShapeOfArray & s,long nn): n(s.n),step(s.n),next(nn) {}
  ShapeOfArray(long nn): n(nn),step(1),next(-1) {}

  ShapeOfArray(long nn,long s): n(nn),step(s),next(-1) {}

  ShapeOfArray(long nn,long s,long nextt): n(nn),step(s),next(nextt) {}

  ShapeOfArray(const ShapeOfArray &old,const SubArray &sub)
         : n(sub.n),step(old.step*sub.step),next(old.next)
          { K_throwassert((sub.last())*old.step <= old.last());} // a constructor

  ShapeOfArray(const ShapeOfArray &old,long stepo,long start)
         : n(old.n-start),step(old.step*stepo),next(old.next)
          { K_throwassert(n>=0);}

  long end()  const      { return n*step;}
  long last()      const      { return (n-1)*step;}
  long constant()  const { return step==0;}
  long index(long k) const { K_throwassert( (k>=0) && ( (k <n) || !step) );
           return step*k;}
  ShapeOfArray operator*(long stepp) const {return  ShapeOfArray(n,step*stepp,next);}
  bool SameShape(const ShapeOfArray & a) const
          { return  !step || !a.step || a.n == n ;}
  long  N(const ShapeOfArray & a) { return step ? n : a.n;} // size of 2 shape


// protected:
  long operator[](long k) const {
    //if( k<0 || ( k<n && !step) )
    //  cout << "k,n,step=" << k << " " << n << " " << step << endl;
    K_throwassert( (k>=0) && ( (k <n) || !step) );
           return step*k;}
  void init(long nn,long s=1,long nextt=-1) { n=nn; step=s; next=nextt;}
};

ostream & operator<<(ostream & f,const ShapeOfArray & s);

inline bool  SameShape(const ShapeOfArray & a,const ShapeOfArray & b)
           { return  !a.step || !b.step || a.n == b.n ;}

inline long N(const ShapeOfArray & a,const ShapeOfArray & b)
           { K_throwassert(SameShape(a,b)); return  a.step ? a.n :  b.n ;}

inline SubArray::SubArray(const ShapeOfArray & s)
           :n(s.n),step(s.step),start(0) {}


template<class R>
ostream & operator<<(ostream & f,const KN_<const_R> & v) ;

template<class R> istream & operator>>(istream & f, KN_<R> & v);
template<class R> istream & operator>>(istream & f, KN<R> & v);

template<class R>
class SetArray { public:
    R o,step;
    long n;
    explicit SetArray(long nn,R oo=R(),R sstep=R(1)): o(oo),n(nn),step(sstep) {}
    template<class K>   SetArray(SetArray<K> sa): o(sa.o),n(sa.n),step(sa.step) {}

  R operator[](long i) const { return i <= n ? o + R(i)*step : R();}
    long size() const {return n;}
};

// add for gmm++  march 2010 ....
  template<typename T> struct KN__iterator {
    typedef T                   value_type;
    typedef value_type*         pointer;
    typedef value_type&         reference;
    typedef size_t              size_type;
    typedef ptrdiff_t           difference_type;
    typedef std::bidirectional_iterator_tag iterator_category;
    typedef KN__iterator<T> iterator;
  private:
    T* v;
    long step;
  public:
    reference operator *() const { return *v; }
    pointer operator->() const { return &(operator*()); }

    iterator &operator ++() { v += step; return *this; }
    iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
    iterator &operator --() { v -= step; return *this; }
    iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }

    bool operator ==(const iterator &i) const { return v == i.v; }
    bool operator !=(const iterator &i) const { return !(v == i.v); }

    KN__iterator(T *vv,int s) : v(vv),step(s) {K_throwassert(v);}
    friend class KN__const_iterator<T>;
  };

template<typename T> struct KN__const_iterator {
  typedef T                   value_type;
  typedef const value_type*   pointer;
  typedef const value_type&   reference;
  typedef size_t              size_type;
  typedef ptrdiff_t           difference_type;
  typedef std::forward_iterator_tag iterator_category;
  typedef KN__const_iterator<T> iterator;

private:
  const T* v;
  long step;
public:

    reference operator *() const { return *v; }
    pointer operator->() const { return &(operator*()); }

    iterator &operator ++() { v += step; return *this; }
    iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
    iterator &operator --() { v -= step; return *this; }
    iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }

    bool operator ==(const iterator &i) const { return v == i.v; }
    bool operator !=(const iterator &i) const { return !(v == i.v); }

  KN__const_iterator(const T *vv,int s) : v(vv),step(s) {K_throwassert(v);}
  KN__const_iterator(const KN__const_iterator<T> & i) : v(i.v),step(i.step){}

};


template<class R>
class KN_: public  ShapeOfArray {
protected:
  R *v;
public:
  typedef R K; // type of data
  typedef KN__iterator<R> iterator;
  typedef KN__const_iterator<R> const_iterator;
  iterator end(){ return iterator(v+n*step,step);}
  iterator begin(){ return iterator(v,step);}
  const_iterator end() const { return const_iterator(v+n*step,step);}
  const_iterator begin() const { return const_iterator(v,step);}
  void clear() { *this=R();}
  long N() const {return n;}
  bool unset() const { return !v;}
  void set(R * vv,int nn,int st=1,int nx=-1) {v=vv;n=nn;step=st;next=nx;}
  long size() const{return step?n*step:n;}
  operator R *() const {return v;}
  KN_(const KN_<R> & u) :ShapeOfArray(u),v(u.v){}
  KN_(const KN_<R> & U,const SubArray & sa)  : ShapeOfArray(U,sa),v(U.v + U.index(sa.start)) {}

  KN_ operator()(const SubArray & sa) const { return KN_(*this,sa);} // sub array

  R & operator[](long i) const {return v[index(i)];}
  R & operator()(long i) const {return v[index(i)];}
  R & operator[](int i) const {return v[index(i)];}
  R & operator()(int i) const {return v[index(i)];}

  R operator,(const KN_<const_R> & v) const; // dot  product

   KN_& operator  =(const SetArray<R> & u)  ;
   KN_& operator +=(const SetArray<R> & u)  ;
   KN_& operator -=(const SetArray<R> & u)  ;
   KN_& operator *=(const SetArray<R> & u)  ;
   KN_& operator /=(const SetArray<R> & u)  ;

    KN_& operator  =(const KN_<const_R> & u)  ;
    KN_& operator +=(const KN_<const_R> & u)  ;
    KN_& operator -=(const KN_<const_R> & u)  ;

    KN_& operator *=(const KN_<const_R> & u)  ;
    KN_& operator /=(const KN_<const_R> & u)  ;


   KN_& operator = (const_R  a) ;
   KN_& operator +=(const_R  a) ;
   KN_& operator -=(const_R  a) ;
   KN_& operator /=(const_R  a) ;
   KN_& operator *=(const_R  a) ;

   KN_& operator  = (R*  a) { return operator =(KN_<R>(a,n));}
   KN_& operator += (R*  a) { return operator+=(KN_<R>(a,n));}
   KN_& operator -= (R*  a) { return operator-=(KN_<R>(a,n));}
   KN_& operator *= (R*  a) { return operator*=(KN_<R>(a,n));}
   KN_& operator /= (R*  a) { return operator/=(KN_<R>(a,n));}


  R min() const ;
  R max() const ;
  R sum() const ;
  double norm() const ;
  double l2() const ;
  double l1() const ;
  double linfty() const ;
  double lp(double p) const ;

  template<class T> long last(const T &) const;
  template<class T> long first(const T &) const;

  void map(R (*f)(R )); // apply the f fonction a all element of the array
  void map(R (*f)(const  R& )); // apply the f fonction a all element of the array

 template<class T>
   void set(R (*f)(const  T& ),KN_<T> & u); // apply the f fonction a all element of the array u

   KN_& operator =(const DotStar_KN_<R> & u) ;
   KN_& operator+=(const DotStar_KN_<R> & u) ;
   KN_& operator-=(const DotStar_KN_<R> & u) ;
   KN_& operator*=(const DotStar_KN_<R> & u) ;
   KN_& operator/=(const DotStar_KN_<R> & u) ;

   KN_& operator =(const DotSlash_KN_<R> & u) ;
   KN_& operator+=(const DotSlash_KN_<R> & u) ;
   KN_& operator-=(const DotSlash_KN_<R> & u) ;
   KN_& operator*=(const DotSlash_KN_<R> & u) ;
   KN_& operator/=(const DotSlash_KN_<R> & u) ;

   KN_& operator =(const if_KN_<R> & u) ;
   KN_& operator+=(const if_KN_<R> & u) ;
   KN_& operator-=(const if_KN_<R> & u) ;
   KN_& operator*=(const if_KN_<R> & u) ;
   KN_& operator/=(const if_KN_<R> & u) ;

   KN_& operator =(const ifnot_KN_<R> & u) ;
   KN_& operator+=(const ifnot_KN_<R> & u) ;
   KN_& operator-=(const ifnot_KN_<R> & u) ;
   KN_& operator*=(const ifnot_KN_<R> & u) ;
   KN_& operator/=(const ifnot_KN_<R> & u) ;

   KN_& operator =(const Add_KN_<R> & u) ;
   KN_& operator+=(const Add_KN_<R> & u) ;
   KN_& operator-=(const Add_KN_<R> & u) ;
   KN_& operator*=(const Add_KN_<R> & u) ;
   KN_& operator/=(const Add_KN_<R> & u) ;

   template<class I,class T> KN_& operator =  (const KN_ITAB<T,I> & u);
   template<class I,class T> KN_& operator +=  (const KN_ITAB<T,I> & u);
   template<class I,class T> KN_& operator -=  (const KN_ITAB<T,I> & u);
   template<class I,class T> KN_& operator *=  (const KN_ITAB<T,I> & u);
   template<class I,class T> KN_& operator /=  (const KN_ITAB<T,I> & u);


   KN_ITAB< KN_<R>,const KN_<int> >  operator()(const KN_<int> &itab) ;
   KN_ITAB< KN_<R>,const KN_<long> >  operator()(const KN_<long> &itab) ;
   KN_ITAB<const KN_<R>,const KN_<int> > operator()(const KN_<int> &itab) const ;
   KN_ITAB<const KN_<R>,const KN_<long> >  operator()(const KN_<long> &itab) const ;


   KN_& operator =(const Sub_KN_<R> & u) ;
   KN_& operator-=(const Sub_KN_<R> & u) ;
   KN_& operator+=(const Sub_KN_<R> & u) ;
   KN_& operator*=(const Sub_KN_<R> & u) ;
   KN_& operator/=(const Sub_KN_<R> & u) ;

   KN_& operator =(const Mulc_KN_<R> & u) ;
   KN_& operator+=(const Mulc_KN_<R> & u) ;
   KN_& operator-=(const Mulc_KN_<R> & u) ;
   KN_& operator*=(const Mulc_KN_<R> & u) ;
   KN_& operator/=(const Mulc_KN_<R> & u) ;

   KN_& operator =(const Add_Mulc_KN_<R> & u) ;
   KN_& operator+=(const Add_Mulc_KN_<R> & u) ;
   KN_& operator-=(const Add_Mulc_KN_<R> & u) ;
   KN_& operator*=(const Add_Mulc_KN_<R> & u) ;
   KN_& operator/=(const Add_Mulc_KN_<R> & u) ;

   KN_& operator =(const if_arth_KN_<R> & u) ;
   KN_& operator+=(const if_arth_KN_<R> & u) ;
   KN_& operator-=(const if_arth_KN_<R> & u) ;
   KN_& operator*=(const if_arth_KN_<R> & u) ;
   KN_& operator/=(const if_arth_KN_<R> & u) ;


   KN_& operator =(const Mul_KNM_KN_<R> & u) ;
   KN_& operator+=(const Mul_KNM_KN_<R> & u) ;
   KN_& operator-=(const Mul_KNM_KN_<R> & u) ;
   KN_& operator*=(const Mul_KNM_KN_<R> & u) ;
   KN_& operator/=(const Mul_KNM_KN_<R> & u) ;

 //  KN_& operator =(const MatriceCreuseMulKN_<R> & ) ;
 //  KN_& operator +=(const MatriceCreuseMulKN_<R> & ) ;
   KN_& operator =(const typename VirtualMatrice<R>::plusAx & Ax)
    {*this=R(); Ax.A->addMatMul(Ax.x,*this);return *this;}
   KN_& operator =(const typename VirtualMatrice<R>::plusAtx & Ax)
    {*this=R(); Ax.A->addMatTransMul(Ax.x,*this);return *this;}
   KN_& operator +=(const typename VirtualMatrice<R>::plusAx & Ax)
    {  Ax.A->addMatMul(Ax.x,*this);return *this;}
   KN_& operator +=(const typename VirtualMatrice<R>::plusAtx & Ax)
    {  Ax.A->addMatTransMul(Ax.x,*this);return *this;}
   KN_& operator =(const typename VirtualMatrice<R>::solveAxeqb & Ab)
    {*this=R(); Ab.A->Solve(*this,Ab.b);return *this;}

  template<class  A,class B,class C> KN_&  operator =  (const F_KN_<A,B,C>  & u) ;
  template<class  A,class B,class C> KN_&  operator +=  (const F_KN_<A,B,C>  & u) ;
  template<class  A,class B,class C> KN_&  operator -=  (const F_KN_<A,B,C>  & u) ;
  template<class  A,class B,class C> KN_&  operator /=  (const F_KN_<A,B,C>  & u) ;
  template<class  A,class B,class C> KN_&  operator *=  (const F_KN_<A,B,C>  & u) ;


//   KN_& operator =(const MatriceCreuseDivKN_<R> &)  ;

 friend   ostream & operator<< <R>(ostream & f,const KN_<const_R> & v)  ;

  KN_(R *u,const ShapeOfArray & s):ShapeOfArray(s),v(u){}
  KN_(R *u,long nn,long s):ShapeOfArray(nn,s),v(u){}
  KN_(R *u,long nn,long s,long nextt):ShapeOfArray(nn,s,nextt),v(u){}
  KN_(R *u,long nn):ShapeOfArray(nn),v(u){}


  TKN_<R>  t() ; // transpose
  const TKN_<R>  t() const ; // transpose
  notKN_<R>  operator!()  ; //  not
  const  notKN_<R>  operator!() const ; // not

  //  operator KN<R> &();
  // operator const KN<R> &() const;

 private:

  KN_&  operator++(){K_throwassert(next>=0);v += next;return *this;} //    ++U
  KN_&  operator--(){K_throwassert(next>=0);v -= next;return *this;} //    --U
  KN_   operator++(int ){K_throwassert(next>=0); KN_ old=*this;v = v +next;return old;} // U++
  KN_   operator--(int ){K_throwassert(next>=0); KN_ old=*this;v = v -next;return old;} // U++

  KN_(const KN_<R> & u,long offset) :ShapeOfArray(u),v(&u[offset]){}
  KN_(const KN_<R> & u,const ShapeOfArray &sh,long startv=0)
         :ShapeOfArray(sh*u.step),v(&u[startv]){}
  KN_(const KN_<R> & u,long nnext,const ShapeOfArray &sh,long startv=0)
         :ShapeOfArray(sh.n,sh.step*u.step,nnext),v(&u[startv]){ }
  // Add for gmm++ compatibityly ...................
  //  iterator end() {return iterator(v,step)


// friend class KN_<R>;
 friend class KNM_<R>;
 friend class KNMK_<R>;
 friend class KN<R>;
 friend class KNM<R>;
 friend class KNMK<R>;

};

template<class R>
class KNM_: public KN_<R> {
  public:
  ShapeOfArray shapei;
  ShapeOfArray shapej;
  public:
  long IsVector1() const  {  return (shapei.n*shapej.n) ==  this->n ;}
  long N() const {return shapei.n;}
  long M() const {return shapej.n;}
  long size() const { return shapei.n*shapej.n;}

  KNM_(R* u,const ShapeOfArray & s,
            const ShapeOfArray & si,
            const ShapeOfArray & sj)
             : KN_<R>(u,s),shapei(si),shapej(sj){}
  KNM_(R* u,long n,long m)
             : KN_<R>(u,ShapeOfArray(n*m)),shapei(n,1,n),shapej(m,n,1){}
  KNM_(R* u,long n,long m,long s)
             : KN_<R>(u,ShapeOfArray(n*m,s)),shapei(n,1,n),shapej(m,n,1){}
  KNM_(KN_<R> u,long n,long m)
             : KN_<R>(u,ShapeOfArray(m*n)),shapei(n,1,n),shapej(m,n,1){ }

  KNM_(const KN_<R> &u,const ShapeOfArray & si,const ShapeOfArray & sj,long offset=0)
             : KN_<R>(&u[offset],si.last()+sj.last()+1,u.step),shapei(si),shapej(sj)
             {K_throwassert( offset>=0 && this->n+ (this->v-(R*)u) <= u.n);}
  KNM_(const KN_<R> &u,const ShapeOfArray & si,const ShapeOfArray & sj,long offset,long nnext)
             : KN_<R>(&u[offset],si.last()+sj.last()+1,u.step,nnext),shapei(si),shapej(sj)
             {K_throwassert( offset>=0 && this->n+ (this->v-(R*)u) <= u.n);}

  KNM_(KNM_<R> U,const SubArray & si,const SubArray & sj)
             :KN_<R>(U,SubArray(U.ij(si.len1(),sj.len1())+1,U.ij(si.start,sj.start))),
                     shapei(U.shapei,si),shapej(U.shapej,sj){}

  KNM_(KNM_<R> U,const SubArray & sa,const SubArray & si,const SubArray & sj)
             :KN_<R>(U,SubArray(sa)),shapei(U.shapei,si),shapej(U.shapej,sj){}

  KNM_(const KNM_<R> & u)
             :KN_<R>(u),shapei(u.shapei),shapej(u.shapej) {}

  KNM_ operator()(const SubArray & sa,const SubArray & sb) const
            { return KNM_(*this,sa,sb);} // sub array

  long ij(long i,long j) const
            { return shapei.index(i)+shapej.index(j);}
  long indexij(long i,long j)        const
            { return this->index(shapei.index(i)+shapej.index(j));}
  R & operator()(long i,long j)     const
            { return this->v[indexij(i,j)];}
  R & operator()(int i,int j)     const
            { return this->v[indexij(i,j)];}


  KN_<R> operator()(const char,long j    )  const   // une colonne j  ('.',j)
            { return KN_<R>(&this->v[this->index(shapej.index(j))],shapei*this->step);}
  KN_<R> operator()(long i    ,const char)  const   // une ligne i  (i,'.')
            { return KN_<R>(&this->v[this->index(shapei.index(i))],shapej*this->step);}
  KN_<R> operator()(const char,int j    )  const   // une colonne j  ('.',j)
            { return KN_<R>(&this->v[this->index(shapej.index(j))],shapei*this->step);}
  KN_<R> operator()(int i    ,const char)  const   // une ligne i  (i,'.')
            { return KN_<R>(&this->v[this->index(shapei.index(i))],shapej*this->step);}
  KN_<R> operator()(const char,const char)  const   // tous
            { return *this;}
  KNM_<R> t() const
    { return KNM_<R>(this->v,*this,shapej,shapei);} //  before  { return KNM_<R>(*this,shapej,shapei,v);}

   KNM_& operator =(const KNM_<const_R> & u) ;
   KNM_& operator =(const_R a)               ;
   KNM_& operator+=(const_R a)               ;
   KNM_& operator-=(const_R a)               ;
   KNM_& operator/=(const_R a)               ;
   KNM_& operator*=(const_R a)               ;
   KNM_& operator+=(const KNM_<const_R> & u) ;
   KNM_& operator-=(const KNM_<const_R> & u) ;
   KNM_& operator*=(const KNM_<const_R> & u) ;
   KNM_& operator/=(const KNM_<const_R> & u) ;

   KNM_ &operator =(const outProduct_KN_<R> &);
   KNM_ &operator +=(const outProduct_KN_<R> &);
   KNM_ &operator -=(const outProduct_KN_<R> &);
   KNM_ &operator /=(const outProduct_KN_<R> &); // bofbof
   KNM_ &operator *=(const outProduct_KN_<R> &); // bofbof

private:
  KNM_& operator++() {this->v += this->next;return *this;} // ++U
  KNM_& operator--() {this->v -= this->next;return *this;} // ++U
  KNM_  operator++(int ){KNM_<R> old=*this;this->v = this->v +this->next;return old;} // U++
  KNM_  operator--(int ){KNM_<R> old=*this;this->v = this->v -this->next;return old;} // U--


 friend class KN_<R>;
// friend class KNM_<R>;
 friend class KNMK_<R>;
 friend class KN<R>;
 friend class KNM<R>;
 friend class KNMK<R>;
};

template<class T,class I>
struct KN_ITAB
{
  KN_ITAB(const T &vv,const I &iindex) : v(vv),index(iindex) {}
  T  v;
  I  index;
  KN_ITAB & operator=(const T & t);
  KN_ITAB & operator+=(const T & t);
  KN_ITAB & operator-=(const T & t);
  KN_ITAB & operator*=(const T & t);
  KN_ITAB & operator/=(const T & t);
  typename T::R & operator[](long i){ return v[index[i]];}
  const typename T::R &  operator[](long i) const { return v[index[i]];}
  long N() const { return index.N();}
};

template<class R>  KN_ITAB<const KN_<R>,const KN_<int> >  KN_<R>::operator()(const KN_<int>  &itab) const { return KN_ITAB<const KN_<R>,const KN_<int> > (*this,itab);}
template<class R>  KN_ITAB<const KN_<R>,const KN_<long> >  KN_<R>::operator()(const KN_<long>  &itab) const { return KN_ITAB<const KN_<R>,const KN_<long> > (*this,itab);}
template<class R>  KN_ITAB< KN_<R>,const KN_<int> >  KN_<R>::operator()(const KN_<int>  &itab) { return KN_ITAB<KN_<R>,const KN_<int> > (*this,itab);}
template<class R>  KN_ITAB< KN_<R>,const KN_<long> >  KN_<R>::operator()(const KN_<long>  &itab) { return KN_ITAB<KN_<R>,const KN_<long> > (*this,itab);}


template<class R>
struct TKN_:public KN_<R> {
    TKN_(const KN_<R> &x) : KN_<R>(x) {}
};

template<class R>
struct notKN_:public KN_<R> {
    notKN_(const KN_<R> &x) : KN_<R>(x) {}
    notnotKN_<R>  operator!()  ; //  not
    const  notnotKN_<R>  operator!() const ; // not
};

template<class R>
struct notnotKN_:public KN_<R> {
    notnotKN_(const notKN_<R> &x) : KN_<R>(x) {}
    notKN_<R>  operator!()  ; //  notnot
    const  notKN_<R>  operator!() const ; // notnot
};

template<class R>
TKN_<R>  KN_<R>::t() { return *this;} // transpose

template<class R>
const TKN_<R>  KN_<R>::t() const { return *this;} // transpose

template<class R>
notKN_<R>  KN_<R>::operator!() { return *this;} // not

template<class R>
const notKN_<R>  KN_<R>::operator!() const { return *this;} // not

template<class R>
notnotKN_<R>  notKN_<R>::operator!() { return *this;} // not

template<class R>
const notnotKN_<R>  notKN_<R>::operator!() const { return *this;} // not


template<class R>
struct outProduct_KN_ {
    const KN_<R>  a,b;
    R c;
    long N() const {return a.N();    }
    long M() const {return b.N();    }
    outProduct_KN_(const KN_<R> & aa, const KN_<R> &bb,R cc=(R)1) : a(aa),b(bb),c(cc) {}
    outProduct_KN_(const KN_<R> * aa, const KN_<R> &bb,R cc=(R)1) : a(*aa),b(bb),c(cc) {}
    outProduct_KN_(const KN_<R> * aa, const KN_<R> *bb,R cc=(R)1) : a(*aa),b(*bb),c(cc) {}
    outProduct_KN_(const Mulc_KN_<R> & aa,const KN_<R> & bb) : a(aa.a),b(bb),c(aa.b) {}
    outProduct_KN_ operator * (R cc) { return outProduct_KN_(a,b,c*cc);}
};

template<class R>
struct if_KN_ {
    const KN_<R> & a,&b;
    R c;
    if_KN_(const KN_<R> & aa, const KN_<R> &bb,R cc=1.) : a(aa),b(bb),c(cc) {}
    if_KN_ operator * (R cc) { return if_KN_(a,b,c*cc);}
};

template<class R>
struct ifnot_KN_ {
    const KN_<R> & a,&b;
    R c;
    ifnot_KN_(const KN_<R> & aa, const KN_<R> &bb,R cc=1.) : a(aa),b(bb),c(cc) {}
    ifnot_KN_ operator * (R cc) { return ifnot_KN_(a,b,c*cc);}
};


template<class R>
outProduct_KN_<R> operator*(const KN_<R> &a,const TKN_<R> &b)
{ return outProduct_KN_<R>(a,b);}

template<class R>
ifnot_KN_<R> operator*(const KN_<R> &a,const notKN_<R> &b)
{ return ifnot_KN_<R>(b,a);}

template<class R>
ifnot_KN_<R> operator*(const KN_<R> &a,const notnotKN_<R> &b)
{ return if_KN_<R>(b,a);}

template<class R>
ifnot_KN_<R> operator*(const notKN_<R> &b,const KN_<R> &a)
{ return ifnot_KN_<R>(b,a);}

template<class R>
ifnot_KN_<R> operator*(const notnotKN_<R> &b,const KN_<R> & a)
{ return if_KN_<R>(b,a);}


template<class R>
R operator*(const TKN_<R> &a,const KN_<R> &b)
{ return (a,b);}

template<class R>
class KNMK_: public KN_<R> {
  friend class KNMK<R>;
  public:
  ShapeOfArray shapei;
  ShapeOfArray shapej;
  ShapeOfArray shapek;
  public:
  long IsVector1() const {  return (shapei.n*shapej.n*shapek.n) == this->n ;}
  long N() const {return shapei.n;}
  long M() const {return shapej.n;}
  long K() const {return shapek.n;}
  long size() const { return shapei.n*shapej.n*shapek.n;}
  KNMK_(const ShapeOfArray & s,
        const ShapeOfArray & si,
        const ShapeOfArray & sj,
        const ShapeOfArray & sk,
	    R * u)
    : KN_<R>(u,s),shapei(si),shapej(sj),shapek(sk){}

  KNMK_(R* u,long n,long m,long k)
    : KN_<R>(u, ShapeOfArray(n*m*k)),shapei(n,1,n),shapej(m,n,1),shapek(k,n*m,n*m){};

//  KNMK_(const KN_<R> & u,long n,long m,long k)
//   : KN_<R>(ShapeOfArray(n*m*k)),shapei(n,1,n),shapekj(m,n,1),u),
//     shapek(k,n*m,n*m){};

  KNMK_(const KNMK_<R> &U,const SubArray & si,const SubArray & sj,const SubArray & sk)  :
    KN_<R>(U,SubArray(U.ijk(si.len1(),sj.len1(),sk.len1())+1,
                       U.ijk(si.start,sj.start,sk.start))),
                       shapei(U.shapei,si),
                       shapej(U.shapej,sj),
                       shapek(U.shapek,sk){}

  KNMK_(const KNMK_<R> & u) :KN_<R>(u),shapei(u.shapei),shapej(u.shapej),shapek(u.shapek) {}


  long ijk(long i,long j,long k) const
              { return shapei.index(i)+shapej.index(j)+shapek.index(k);}
  long indexijk(long i,long j,long k) const
              {return this->index(shapei.index(i)+shapej.index(j)+shapek.index(k));}

  R & operator()(long i,long j,long k)   const   {return this->v[indexijk(i,j,k)];}
  R & operator()(int i,int j,int k)   const   {return this->v[indexijk(i,j,k)];}

//  pas de tableau suivant
 KN_<R>  operator()(const char ,long j,long k)  const  { // le tableau (.,j,k)
        return KN_<R>(*this,-1,shapei,shapej[j]+shapek[k]);}
 KN_<R>  operator()(long i,const char ,long k)  const  { // le tableau (i,.,k)
        return KN_<R>(*this,-1,shapej,shapei[i]+shapek[k]);}
 KN_<R>  operator()(long i,long j,const char )  const  { // le tableau (i,j,.)
        return KN_<R>(*this,-1,shapek,shapei[i]+shapej[j]);}

 KN_<R>  operator()(const char ,int j,int k)  const  { // le tableau (.,j,k)
        return KN_<R>(*this,-1,shapei,shapej[j]+shapek[k]);}
 KN_<R>  operator()(int i,const char ,int k)  const  { // le tableau (i,.,k)
        return KN_<R>(*this,-1,shapej,shapei[i]+shapek[k]);}
 KN_<R>  operator()(int i,int j,const char )  const  { // le tableau (i,j,.)
        return KN_<R>(*this,-1,shapek,shapei[i]+shapej[j]);}
//
 KNM_<R>  operator()(const char ,const char ,long k)  const  { // le tableau (.,.,k)
        return KNM_<R>(*this,shapei,shapej,shapek[k],shapek.next);} // step = n*m
 //attention les suivants ne marche pas
 KNM_<R>  operator()(const char ,long j,const char )  const  { // le tableau (.,j,.)
        return KNM_<R>(*this,shapei,shapek,shapej[j],-1/*shapej.next*/);} // step = n

 KNM_<R>  operator()(long i,const char ,const char )  const  { // le tableau (i,.,.)
        return KNM_<R>(*this,shapej,shapek,shapei[i],-1/*shapei.next*/);}  // step = 1

 KNM_<R>  operator()(const char ,const char ,int k)  const  { // le tableau (.,.,k)
        return KNM_<R>(*this,shapei,shapej,shapek[k],shapek.next);} // step = n*m
 //attention les suivants ne marche pas
 KNM_<R>  operator()(const char ,int j,const char )  const  { // le tableau (.,j,.)
        return KNM_<R>(*this,shapei,shapek,shapej[j],-1/*shapej.next*/);} // step = n

 KNM_<R>  operator()(int i,const char ,const char )  const  { // le tableau (i,.,.)
        return KNM_<R>(*this,shapej,shapek,shapei[i],-1/*shapei.next*/);}  // step = 1

   KNMK_& operator =(const KNMK_<const_R> & u) ;
   KNMK_& operator+=(const KNMK_<const_R> & u)  ;
   KNMK_& operator-=(const KNMK_<const_R> & u)  ;
   KNMK_& operator/=(const KNMK_<const_R> & u)  ;
   KNMK_& operator*=(const KNMK_<const_R> & u)  ;
   KNMK_& operator =(const_R a)  ;
   KNMK_& operator+=(const_R a)  ;
   KNMK_& operator-=(const_R a)  ;
   KNMK_& operator/=(const_R a)  ;
   KNMK_& operator*=(const_R a)  ;

  KNMK_  operator()(SubArray si,SubArray sj,SubArray sk) const
        {return KNMK_(*this,si,sj,sk);}

  private:
//  KNMK_&  operator++(){v += next;return *this;} // ++U
//  KNMK_&  operator--(){v -= next;return *this;} // --U
//  KNMK_  operator++(long ){KNMK_ old=*this;v = v +next;return old;} // U++
//  KNMK_  operator--(long ){KNMK_ old=*this;v = v -next;return old;} // U--


friend class KNM_<R>;
friend class KN_<R>;

};


template<class R>
class KN :public KN_<R> { public:

  typedef R K;

 // explicit  KN(const R & u):KN_<R>(new R(uu),1,0) {}
  KN() : KN_<R>(0,0) {}
  KN(long nn) : KN_<R>(new R[nn],nn)         {}
  KN(long nn, R * p) : KN_<R>(new R[nn],nn)
    { KN_<R>::operator=(KN_<R>(p,nn));}
  KN(long nn,R (*f)(long i) ) : KN_<R>(new R[nn],nn)
        {for(long i=0;i<this->n;i++) this->v[i]=f(i);}
  KN(long nn,const  R & a) : KN_<R>(new R[nn],nn)
        { KN_<R>::operator=(a);}
  KN(long nn,long s,const  R  a) : KN_<R>(new R[nn],nn,s)
        { KN_<R>::operator=(a);}
  template<class S>   KN(const KN_<S> & s):KN_<R>(new R[s.n],s.n)
        {for (long i=0;i<this->n;i++) this->v[i] = s[i];}
  template<class S>  KN(const KN_<S> & s,R (*f)(S )):KN_<R>(new R[s.n],s.n)
        {for (long i=0;i<this->n;i++) this->v[i] = f(s[i]);}
  KN(const KN<R> & u):KN_<R>(new R[u.n],u.n)
        { KN_<R>::operator=(u);}
  KN(bool ,KN<R> & u):KN_<R>(u) {u.v=0;u.n=0;}// remove copy for return of local KN.

  //  explicit KN(const KN_<R> & u):KN_<R>(new R[u.n],u.n)
  //      { KN_<R>::operator=(u);}

  ~KN(){delete [] this->v;}

  void CheckSet() { if(!(this->n)) {cerr << "Error RNM set array\n";K_throwassert(0); exit(1);}}
   KN& operator  = (R*  a) { CheckSet(); return operator =(KN_<R>(a,this->n));}
   KN& operator += (R*  a) { CheckSet(); return operator+=(KN_<R>(a,this->n));}
   KN& operator -= (R*  a) { CheckSet(); return operator-=(KN_<R>(a,this->n));}
   KN& operator *= (R*  a) { CheckSet(); return operator*=(KN_<R>(a,this->n));}
   KN& operator /= (R*  a) { CheckSet(); return operator/=(KN_<R>(a,this->n));}

   KN& operator  =(const SetArray<R> & u)
     { if(this->unset())
	 this->set(new R[u.size()],u.size(),0,0);
       KN_<R>::operator= (u);
       return *this;}
   KN& operator +=(const SetArray<R> & u)
     { if(this->unset())
	 this->set(new R[u.size()],u.size(),0,0);
       KN_<R>::operator+= (u);
       return *this;}
   KN& operator -=(const SetArray<R> & u)
     { if(this->unset())  this->set(new R[u.size()],u.size(),0,0); KN_<R>::operator-= (u);return *this;}
   KN& operator *=(const SetArray<R> & u)
     { if(this->unset())  this->set(new R[u.size()],u.size(),0,0); KN_<R>::operator*= (u);return *this;}
   KN& operator /=(const SetArray<R> & u)
     { if(this->unset())  this->set(new R[u.size()],u.size(),0,0); KN_<R>::operator/= (u);return *this;}

   KN& operator =(const_R a)
        { if(this->unset()) this->set(new R[1],1,0,0); KN_<R>::operator= (a);return *this;}
   KN& operator =(const KN_<R>& a)
        { if(this->unset())  this->set(new R[a.N()],a.N()); KN_<R>::operator= (a);return *this;}
   KN& operator =(const KN<R>& a)
        { if(this->unset()) this->set(new R[a.N()],a.N()); KN_<R>::operator= (a);return *this;}
   KN& operator =(const Add_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const DotStar_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const if_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const ifnot_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const DotSlash_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const Sub_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const Mulc_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const Add_Mulc_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}
   KN& operator =(const if_arth_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.a.N()],u.a.N());KN_<R>::operator=(u);return *this;}


   KN& operator =(const Mul_KNM_KN_<R> & u)
        { if(this->unset())  this->set(new R[u.b.N()],u.b.N());KN_<R>::operator=(u);return *this;}
//   KN& operator =(const MatriceCreuseMulKN_<R> & Ax)
//       {if(this->unset())  this->set(new R[Ax.v.N()],Ax.v.N()); KN_<R>::operator=(Ax);return *this;}
//   KN& operator +=(const MatriceCreuseMulKN_<R> & Ax)
//       {if(this->unset())  this->set(new R[Ax.v.N()],Ax.v.N()); KN_<R>::operator+=(Ax);return *this;}
//   KN& operator =(const MatriceCreuseDivKN_<R> & A1x)
//       { if(this->unset())  this->set(new R[A1x.v.N()],A1x.v.N());KN_<R>::operator=(A1x);return *this;}
  // correcton aout 2007 FH  add N,M flied in VirtualMatrice
   KN& operator =(const typename VirtualMatrice<R>::plusAx & Ax)
        { if(this->unset() && Ax.A->N )  this->set(new R[Ax.A->N],Ax.A->N);KN_<R>::operator=(Ax);return *this;}
   KN& operator =(const typename VirtualMatrice<R>::solveAxeqb & Ab)
        { if(this->unset())  this->set(new R[Ab.b.N()],Ab.b.N());KN_<R>::operator=(Ab);return *this;}
   KN& operator +=(const typename  VirtualMatrice<R>::plusAx & Ax)
  { if(this->unset()  && Ax.A->N) {
         this->set(new R[Ax.A->N],Ax.A->N);
        KN_<R>::operator=(R());}
    KN_<R>::operator+=(Ax);
    return *this;}
   KN& operator =(const typename VirtualMatrice<R>::plusAtx & Ax)
        { if(this->unset()&&Ax.A->M)  this->set(new R[Ax.A->M],Ax.A->M);KN_<R>::operator=(Ax);return *this;}
   KN& operator +=(const typename VirtualMatrice<R>::plusAtx & Ax)
  { if(this->unset()&&Ax.A->M) {
        this->set(new R[Ax.A->M],Ax.A->M);
      KN_<R>::operator=(R());}
      KN_<R>::operator+=(Ax);
     return *this;}
// end correcton FH
   template<class P,class Q>
     KN& operator =(const  PplusQ<P,Q> & PQ)
      { *this=PQ.p; *this+=PQ.q;return *this; }
   template<class P,class Q>
     KN& operator +=(const  PplusQ<P,Q> & PQ)
      { *this+=PQ.p; *this+=PQ.q;return *this; }

   KN& operator -=(const_R a)
        { KN_<R>::operator-=(a);return *this;}
   KN& operator -=(const KN_<R>& a)
        { KN_<R>::operator-= (a);return *this;}
   KN& operator -=(const Add_KN_<R> & u)
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const DotStar_KN_<R> & u)
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const DotSlash_KN_<R> & u)
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const Sub_KN_<R> & u)
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const Mulc_KN_<R> & u)
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const Add_Mulc_KN_<R> & u)
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const if_arth_KN_<R> & u)
        { KN_<R>::operator-=(u);return *this;}
   KN& operator -=(const Mul_KNM_KN_<R> & u)
        { KN_<R>::operator-=(u);return *this;}

   KN& operator +=(const_R a)
        { KN_<R>::operator += (a);return *this;}
   KN& operator += (const KN_<R>& a)
        { KN_<R>::operator+= (a);return *this;}
   KN& operator +=(const Add_KN_<R> & u)
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const DotStar_KN_<R> & u)
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const DotSlash_KN_<R> & u)
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const Sub_KN_<R> & u)
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const Mulc_KN_<R> & u)
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const Add_Mulc_KN_<R> & u)
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const if_arth_KN_<R> & u)
        { KN_<R>::operator+=(u);return *this;}
   KN& operator +=(const Mul_KNM_KN_<R> & u)
        { KN_<R>::operator+=(u);return *this;}


   KN& operator/=(const_R a)
        { KN_<R>::operator/=(a);return *this;}
   KN& operator /= (const KN_<R>& a)
        { KN_<R>::operator/= (a);return *this;}
   KN& operator /=(const Add_KN_<R> & u)
        { KN_<R>::operator/=(u);return *this;}
   KN& operator /=(const Sub_KN_<R> & u)
        { KN_<R>::operator/=(u);return *this;}
   KN& operator /=(const Mulc_KN_<R> & u)
        { KN_<R>::operator/=(u);return *this;}
   KN& operator /=(const Add_Mulc_KN_<R> & u)
        { KN_<R>::operator/=(u);return *this;}
   KN& operator /=(const if_arth_KN_<R> & u)
        { KN_<R>::operator/=(u);return *this;}

   KN& operator /=(const Mul_KNM_KN_<R> & u)
        { KN_<R>::operator/=(u);return *this;}

   KN& operator*=(const_R a)
        { KN_<R>::operator*=(a);return *this;}
   KN& operator*=(const KN_<const_R>& a)
        { KN_<R>::operator*= (a);return *this;}
   KN& operator *=(const Add_KN_<R> & u)
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const Sub_KN_<R> & u)
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const Mulc_KN_<R> & u)
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const Add_Mulc_KN_<R> & u)
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const if_arth_KN_<R> & u)
        { KN_<R>::operator*=(u);return *this;}
   KN& operator *=(const Mul_KNM_KN_<R> & u)
        { KN_<R>::operator*=(u);return *this;}


  template<class I,class T> KN& operator =   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator =(ui);  return *this;}
  template<class I,class T> KN& operator +=   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator +=(ui);  return *this;}
  template<class I,class T> KN& operator -=   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator -=(ui);  return *this;}
  template<class I,class T> KN& operator *=   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator *=(ui);  return *this;}
  template<class I,class T> KN& operator /=   (const KN_ITAB<T ,I> & ui)
     {  KN_<R>::operator /=(ui);  return *this;}


  //  two opertor to cast to an array of constant
//    operator KN_<const_R> & ()
//          { return *  (KN_<const_R>*) this;}
//    operator KN_<const_R> const & ()  const
//          { return *(const KN_<const_R>*) this;}
//    operator KN<const_R> & ()
//          { return   (KN<const_R> &) *this;}
//    operator KN<const_R> const & ()  const
//          { return (const KN<const_R>& ) *this;}
  void init(long nn) {this->n=nn;this->step=1;this->next=-1;this->v=new R[nn];}
  void init() {this->n=0;this->step=1;this->next=-1;this->v=0;}
  void init(const KN_<R> & a){init(a.N()); operator=(a);}
  void resize(long nn) {
    if ( nn != this->n)
     {
       R *vo=this->v;
       long no=std::min(this->n,nn), so=this->step;
       ShapeOfArray::init(nn);
       this->v=new R[this->n];
       // copy
       if(this->v && vo)
         for(long i=0,j=0;j<no;i++,j+=so)
           this->v[i]=vo[j];
        delete [] vo;} }//  mars 2010
    void destroy(){ if(this->next++ ==-1) {delete [] this->v; this->v=0;this->n=0;}}//  mars 2010
    void increment() { this->next--;}
};

//  Array with 2 indices
//  ---------------------

template<class R>
class KNM: public KNM_<R>{ public:

  KNM(long nn,long mm)
        :KNM_<R>(new R[nn*mm],nn,mm){}
   KNM(const KNM<R> & u)  // PB si stepi ou stepj nulle
        :KNM_<R>(new R[u.size()],u.N(),u.M())
       { KN_<R>::operator=(u);}
  explicit KNM(const KNM_<R> & u)
        :KNM_<R>(new R[u.size()],u.N(),u.M())
        { KNM_<R>::operator=(u);}

  ~KNM(){delete [] this->v;}

   KNM& operator=(const KNM_<const_R> & u)
    { if(this->unset()) this->init(u.N(),u.M()) ; KNM_<R>::operator=(u);return *this;}
   KNM& operator=(const_R a)
    { if(this->unset()) RNM_FATAL_ERROR(" KNM operator=(double)"); KNM_<R>::operator=(a);return *this;}
   KNM& operator+=(const_R  a)
        { if(this->unset()) RNM_FATAL_ERROR(" KNM operator+=(double)"); KNM_<R>::operator+=(a);return *this;}
   KNM& operator-=(const_R  a)
        {if(this->unset()) RNM_FATAL_ERROR(" KNM operator-=(double)"); KNM_<R>::operator-=(a);return *this;}
   KNM& operator/=(const_R  a)
        {if(this->unset()) RNM_FATAL_ERROR(" KNM operator/=(double)"); KNM_<R>::operator/=(a);return *this;}
   KNM& operator*=(const_R  a)
        {if(this->unset()) RNM_FATAL_ERROR(" KNM operator*=(double)"); KNM_<R>::operator*=(a);return *this;}
   KNM& operator+=(const KNM_<const_R> & u)
        { if(this->unset()) this->init(u.N(),u.M()) ; KNM_<R>::operator+=(u);return *this;}
   KNM& operator-=(const KNM_<const_R> & u)
        {  if(this->unset()) this->init(u.N(),u.M()) ;KNM_<R>::operator-=(u);return *this;}

   KNM& operator/=(const KNM_<const_R> & u)
        {  if(this->unset()) this->init(u.N(),u.M()) ;KNM_<R>::operator/=(u);return *this;}
   KNM& operator*=(const KNM_<const_R> & u)
        { if(this->unset()) this->init(u.N(),u.M()) ; KNM_<R>::operator*=(u);return *this;}


   KNM &operator  =(const outProduct_KN_<R> & u)
        { if(this->unset()) this->init(u.N(),u.M()) ; KNM_<R>::operator =(u);return *this;}
   KNM &operator +=(const outProduct_KN_<R> & u)
        { if(this->unset()) this->init(u.N(),u.M()) ; KNM_<R>::operator+=(u);return *this;}
   KNM &operator -=(const outProduct_KN_<R> & u)
        { if(this->unset()) this->init(u.N(),u.M()) ; KNM_<R>::operator-=(u);return *this;}
   KNM &operator /=(const outProduct_KN_<R> & u)
        {  if(this->unset()) this->init(u.N(),u.M()) ;KNM_<R>::operator/=(u);return *this;}
   KNM &operator *=(const outProduct_KN_<R> & u)
        {  if(this->unset()) this->init(u.N(),u.M()) ;KNM_<R>::operator*=(u);return *this;}


  //  two opertors to cast to un array of constant
//    operator KNM_<const_R> & ()
//          { return *  (KNM_<const_R>*) this;}
//    operator KNM_<const_R> const & ()  const
//          { return *(const KNM_<const_R>*) this;}

//    operator KNM<const_R> & ()
//          { return *  (KNM<const_R>*) this;}
//    operator KNM<const_R> const & ()  const
//          { return *(const KNM<const_R>*) this;}

    void init() { //  add mars 2010 ...
	this->n=0;this->step=1;this->next=-1;this->v=0;
	this->shapei.init(0);
	this->shapej.init(0);}

  void init(long nn,long mm) {
    ShapeOfArray::init(nn*mm);
    this->shapei.init(nn,1,nn);
    this->shapej.init(mm,nn,1),
    this->v=new R[nn*mm];}

  void resize(long nn,long mm) {
    long kk=nn*mm;

    long lso = this->size();
    long n = this->shapei.n;
    long m = this->shapej.n;

    if( n !=nn && m != mm)
     {
    KNM_ <R> old(*this);
    long no=std::min(n,nn);
    long mo=std::min(m,mm);
    R *vo=this->v;

    // new mat
    ShapeOfArray::init(kk);
    this->v=new R[this->n];
    this->shapei.init(nn,1,nn);
    this->shapej.init(mm,nn,1);

    if(this->v && vo)  // copy
        (*this)(SubArray(no),SubArray(mo)) = old(SubArray(no),SubArray(mo));

    delete []vo;
    }

   }
    void destroy(){ if(this->next++ ==-1) {delete [] this->v; this->v=0;this->n=0;}}
    void increment() {  this->next--;}

//  void destroy(){delete [] this->v;this->n=0 ;}

};

//  Array with 3 indices
//  ---------------------
template<class R>
class KNMK: public KNMK_<R>{ public:

  KNMK(long n,long m,long k)
     :KNMK_<R>(new R[n*m*k],n,m,k){}
  explicit KNMK(const KNMK_<R> & u)
     :KNMK_<R>(new R[u.size()],u.N(),u.M(),u.K())
     { KNMK_<R>::operator=(u);}
   KNMK(const KNMK<R> & u)
     :KNMK_<R>(new R[u.size()],u.N(),u.M(),u.K())
     { KNMK_<R>::operator=(u);}

  ~KNMK(){delete [] this->v;}

  KNMK& operator=(const KNMK_<const_R> & u)
     { KNMK_<R>::operator=(u);return *this;}
  KNMK& operator=(const_R a)
     { KNMK_<R>::operator=(a);return *this;}
  KNMK& operator+=(const_R  a)
     { KNMK_<R>::operator+=(a);return *this;}
  KNMK& operator-=(const_R  a)
     { KNMK_<R>::operator-=(a);return *this;}
  KNMK& operator/=(const_R  a)
     { KNMK_<R>::operator/=(a);return *this;}
  KNMK& operator*=(const_R  a)
     { KNMK_<R>::operator*=(a);return *this;}
  KNMK& operator+=(const KNMK_<const_R> & u)
     { KNMK_<R>::operator+=(u);return *this;}
  // ici jd
  KNMK& operator-=(const KNMK_<const_R> & u)
     { KNMK_<R>::operator-=(u);return *this;}
  KNMK& operator*=(const KNMK_<const_R> & u)
     { KNMK_<R>::operator*=(u);return *this;}
  KNMK& operator/=(const KNMK_<const_R> & u)
     { KNMK_<R>::operator/=(u);return *this;}

//  two opertor to cast to un array of constant
//    operator KNMK_<const_R> & ()
//       { return *  (KNMK_<const_R>*) this;}
//    operator KNMK_<const_R> const & ()  const
//       { return *(const KNMK_<const_R>*) this;}

//    operator KNMK<const_R> & ()
//       { return *  (KNMK<const_R>*) this;}
//    operator KNMK<const_R> const & ()  const
//       { return *(const KNMK<const_R>*) this;}
};

//  -------------  optimization ---------------------
template<class R>
class conj_KN_{public:
  const KN_<const_R> & a;
  conj_KN_(const KN_<const_R> & aa) : a(aa){}
};


inline const KN_<long> conj(const KN_<long> &a){ return a;}
inline const KN_<double> conj(const KN_<double> &a){ return a;}
inline const KN_<float> conj(const KN_<float> &a){ return a;}

//template<class R> conj_KN_<R> conj(const KN<R> &a){ return a;}
template<class R> conj_KN_<R> conj(const KN_<R> &a){ return a;}

template<class R>
class DotStar_KN_{public:
  const KN_<const_R>  a; const KN_<const_R>  b;
  DotStar_KN_(const KN_<const_R> & aa,const KN_<const_R> & bb) : a(aa),b(bb)  {}
 };


template<class R>
class DotSlash_KN_{public:
  const KN_<const_R>  a; const KN_<const_R>  b;
  DotSlash_KN_(const KN_<const_R> & aa,const KN_<const_R> & bb) : a(aa),b(bb)  {}
 };

template<class R>
class Add_KN_{public:
  const KN_<const_R>  a; const KN_<const_R>  b;
  Add_KN_(const KN_<const_R> & aa,const KN_<const_R> & bb)
     : a(aa),b(bb)  { K_throwassert(SameShape(a,b));}
 };

template<class R>
class Sub_KN_{public:
  const KN_<const_R>  a; const KN_<const_R>  b;
  Sub_KN_(const KN_<const_R> & aa,const KN_<const_R> & bb)
    : a(aa),b(bb) { K_throwassert(SameShape(a,b));}
 };

template<class R>
class Mulc_KN_ { public:
  const KN_<const_R>  a;  const_R  b;
  Mulc_KN_(const KN_<const_R> & aa,const_R  bb) : a(aa),b(bb) {}
  Mulc_KN_(const Mulc_KN_<R> & aa,const_R  bb) : a(aa.a),b(aa.b*bb) {}
  Mulc_KN_ operator-() const {return Mulc_KN_(a,-b);}
  outProduct_KN_<R> operator*(const  TKN_<double> & bb)
{  return outProduct_KN_<R>(a,bb,b);}

 };

template<class R>
class Add_Mulc_KN_ { public:
  const KN_<const_R> a,b;
  const R ca,cb;
  Add_Mulc_KN_(const Mulc_KN_<R> & aa,const Mulc_KN_<R> & bb)
        : a(aa.a),b(bb.a),ca(aa.b),cb(bb.b) { K_throwassert(SameShape(a,b));}
  Add_Mulc_KN_(const Mulc_KN_<R> & aa,const KN_<const_R> & bb,const R cbb)
        : a(aa.a),b(bb),ca(aa.b),cb(cbb)  { K_throwassert(SameShape(a,b));}
  Add_Mulc_KN_(const KN_<const_R> & aa,const R caa,const KN_<const_R> & bb,const R cbb)
        : a(aa),b(bb),ca(caa),cb(cbb) { K_throwassert(SameShape(a,b));}
 };

template<class R>
class if_arth_KN_ { public:
  const KN_<const_R> a,b,c;
  if_arth_KN_(const KN_<R> & aa,const KN_<R> & bb,const KN_<R> & cc)
        : a(aa),b(bb),c(cc){ K_throwassert(SameShape(a,b)&&SameShape(a,c));}
 };


template<class R>
class Mul_KNM_KN_ { public:
  const KNM_<const_R> &A;
  const KN_<const_R> &b;
  Mul_KNM_KN_(const  KNM_<const_R>  &aa,const KN_<const_R>  &bb)
        : A(aa),b(bb) {K_throwassert(SameShape(A.shapej,b));}
};


ostream & operator<<(ostream & f,const ShapeOfArray & s);

template<class R> ostream & operator<<(ostream & f,const KN_<const_R>   & v);
template<class R> ostream & operator<<(ostream & f,const KNM_<const_R>  & v);
template<class R> ostream & operator<<(ostream & f,const KNMK_<const_R> & v);
template<class R> inline ostream & operator<<(ostream & f,const KN<const_R>   & v)
    { return f << (const KN_<const_R> &) v;}
template<class R> inline ostream & operator<<(ostream & f,const KNM<const_R>  & v)
    { return f << (const KNM_<const_R> &) v;}
template<class R> inline ostream & operator<<(ostream & f,const KNMK<const_R> & v)
    { return f << (const KNMK_<const_R> &) v;}


template<class R> inline Add_KN_<R> operator+(const KN_<const_R> &a,const KN_<const_R> &b)
    { return Add_KN_<R>(a,b);}
template<class R> inline Sub_KN_<R> operator-(const KN_<const_R> &a,const KN_<const_R> &b)
    { return Sub_KN_<R>(a,b);}
template<class R> inline Mulc_KN_<R> operator*(const KN_<const_R> &a,const R &b)
    { return Mulc_KN_<R>(a,b);}
template<class R> inline Mulc_KN_<R> operator/(const KN_<const_R> &a,const R &b)
    { return Mulc_KN_<R>(a,1/b);}
template<class R> inline Mulc_KN_<R> operator*(const R &b,const KN_<const_R> &a)
    { return Mulc_KN_<R>(a,b);}
template<class R> inline Mulc_KN_<R> operator-(const KN_<const_R> &a)
    { return Mulc_KN_<R>(a,-1);}


template<class R> inline Add_Mulc_KN_<R> operator+(const  Mulc_KN_<R>& a,const Mulc_KN_<R> &b)
    { return Add_Mulc_KN_<R>(a,b);}
template<class R> inline Add_Mulc_KN_<R> operator-(const  Mulc_KN_<R>& a,const Mulc_KN_<R> &b)
    { return Add_Mulc_KN_<R>(a,b.a,-b.b);}

template<class R> inline Add_Mulc_KN_<R> operator+(const  Mulc_KN_<R>& a,const KN_<const_R> &b)
   { return Add_Mulc_KN_<R>(a,b,R(1));}
template<class R> inline Add_Mulc_KN_<R> operator-(const  Mulc_KN_<R>& a,const KN_<const_R> &b)
   { return Add_Mulc_KN_<R>(a,b,R(-1));}

template<class R> inline Add_Mulc_KN_<R> operator+(const KN_<const_R> & b,const  Mulc_KN_<R>& a)
   { return Add_Mulc_KN_<R>(a,b,R(1));}

// modif FH mars 2007
template<class R> inline Add_Mulc_KN_<R> operator-(const KN_<const_R> & a,const  Mulc_KN_<R>& b)
   { return Add_Mulc_KN_<R>(a,R(1),b.a,-b.b);}// modif FH mars 2007

template<class R> inline Mul_KNM_KN_<R> operator*(const  KNM_<const_R> & A,const  KN_<const_R> & b)
    { return Mul_KNM_KN_<R>(A,b);}


template<class R> inline bool  SameShape(const ShapeOfArray & a,const Add_Mulc_KN_<R> & b)
           { return SameShape(a,b.a) ;}
template<class R> inline bool  SameShape(const ShapeOfArray & a,const if_arth_KN_<R> & b)
           { return SameShape(a,b.a) ;}
template<class R> inline bool  SameShape(const ShapeOfArray & a,const Add_KN_<R> & b)
           { return SameShape(a,b.a) ;}
template<class R> inline bool  SameShape(const ShapeOfArray & a,const Sub_KN_<R> & b)
           { return SameShape(a,b.a) ;}
template<class R> inline bool  SameShape(const ShapeOfArray & a,const Mulc_KN_<R> & b)
           { return SameShape(a,b.a) ;}
template<class R> inline bool  SameShape(const ShapeOfArray & a,const DotStar_KN_<R> & b)
           { return SameShape(a,b.a) ;}
template<class R> inline bool  SameShape(const ShapeOfArray & a,const DotSlash_KN_<R> & b)
           { return SameShape(a,b.a) ;}
template<class R> inline bool  SameShape(const ShapeOfArray & a,const Mul_KNM_KN_<R> & b)
           { return a.n==b.A.N() ;}
 inline bool  SameShape(const ShapeOfArray & ,const VirtualMatrice<double>::plusAx & )
           { return true ;} //  pas de test car la matrice peut etre rectangulaire
 inline bool  SameShape(const ShapeOfArray & ,const VirtualMatrice<double>::plusAtx & )
           { return true ;} //  pas de test car la matrice peut etre rectangulaire
 inline bool  SameShape(const ShapeOfArray & ,const VirtualMatrice<complex<double> >::plusAx & )
           { return true ;} //  pas de test car la matrice peut etre rectangulaire
 inline bool  SameShape(const ShapeOfArray & ,const VirtualMatrice<complex<double> >::plusAtx & )
           { return true ;} //  pas de test car la matrice peut etre rectangulaire

 inline bool  SameShape(const ShapeOfArray & ,const double)
           { return true;}
 inline bool  SameShape(const ShapeOfArray & ,const complex<double>)
           { return true;}
 inline bool  SameShape(const ShapeOfArray & ,const complex<float>)
           { return true;}

template<class R>
 inline bool SameShape(KNM<R>& m, const outProduct_KN_<R>& p)
 { return p.a.N()>=m.N() && m.M()>=p.b.N(); }

template<class R> inline long SameAdress(const KN_<R> &a, const KN_<R> &b) { return &a[0]==&b[0];}
// bof -bof
//template<class R> inline
//  KN_<R>::operator KN<R> &() { return *(KN<R> *) (void *) this;}
//template<class R> inline
//  KN_<R>::operator const KN<R> &() const { return *(const KN<R> *) ( const void *) this;}

//  operateur y=Ax-b ou y=Ax + b pour le GC
template<class R>
   PplusQ< typename VirtualMatrice<R>::plusAx, Mulc_KN_<R> >  operator-(const typename VirtualMatrice<R>::plusAx & A,const KN_<R> & B)
    { return PplusQ< typename VirtualMatrice<R>::plusAx, Mulc_KN_<R> >(A,Mulc_KN_<R>(B,R(-1.)));}

template<class R>
   PplusQ< typename VirtualMatrice<R>::plusAx, KN_<R> >  operator+(const typename VirtualMatrice<R>::plusAx & A,const KN_<R> & B)
    { return PplusQ< typename VirtualMatrice<R>::plusAx, KN_<R> >(A,B);}

template<class R>
   PplusQ< typename VirtualMatrice<R>::plusAx, Mulc_KN_<R> >  operator-(const typename VirtualMatrice<R>::plusAx & A,const KN<R> & B)
    { return PplusQ< typename VirtualMatrice<R>::plusAx, Mulc_KN_<R> >(A,Mulc_KN_<R>(B,R(-1.)));}

template<class R>
   PplusQ< typename VirtualMatrice<R>::plusAx, KN_<R> >  operator+(const typename VirtualMatrice<R>::plusAx & A,const KN<R> & B)
    { return PplusQ< typename VirtualMatrice<R>::plusAx, KN_<R> >(A,B);}


template<class R>
KN_<R> diagonal(const KNM<R> & A) {
  K_throwassert(A.N() == A.M());
  return KN_<R>(A,SubArray(A.N(),0,A.N()+1));}

// to def  inv permutation FH mars 2006
class Inv_KN_long{ public:
  KN_<long>  t;
  Inv_KN_long(const KN_<long> & v)
   : t(v) {}
  Inv_KN_long( KN_<long> const  * & v)
   : t(*v) {}
  operator const KN_<long> & () const {return t;}
};


template<class R,typename A,typename B=R> class  F_KN_
{
  public:
  A (*f)(B);
  KN_<R> a;
  long N() const {return a.N();}
  F_KN_( A (*ff)(B),const KN_<R> & aa): f(ff),a(aa) {}
  A operator[](long i) const { return f(a[i]);}
  bool check(long n)  const { return  n <= a.N() || a.constant(); }
  bool constant() const {return a.constant();}
};

template<class R,typename A,typename B>
inline bool  SameShape(const ShapeOfArray & a,const F_KN_<R,A,B>  & b)
           { return  !a.step || b.constant()  || a.n == b.N() ;}

#include "RNM_tpl.hpp"
#ifdef K_throwassert
#undef K_throwassert
#endif
#endif