Newsgroup: it.comp.lang.c++
Titolo: ancora matrici e classi derivate... spero di non rompere troppo

In riferimento alla scrittura di una classe matrice
se definisco la mia classe matrice come
matrix : public valarray
{
int dim1_;
int dim2_;
}
questa eredita le funzioni di valarray ... tra cui per esempio 
l'operatore +.
Questo pero' mi restituisce (giustamente) un valarray...
esiste una maniera molto sintetica ed efficiente di ridefinire 
l'operatore +  sfruttando l'operatore + del valarray già costruito molto 
efficiente?
In pratica la funzione dovrebbe utilizzare l'operatore + del valarray e 
poi aggiungere le dimensioni della matrice dim1_ e dim2_
-- 
Hai recensioni nel cassetto di musica, cinema, letteratura,fumetti??
Ascolta e pentiti e se hai tempo vergognati
http://tonysuper.altervista.org

 poi aggiungere le dimensioni della matrice dim1_ e dim2_
Credo che l'unica sia sfruttare un pattern proxy, e non ci sia nessuna 
trovata piu' furba :/
-- 
Mat/tia Belletti - Undergraduate student @ 
cs.unibo.it
ICQ: 33292311    - email: 
mbel...@cs.unibo.it
IRC: RedGlow     - site(s): 
http://www.cs.unibo.it/~mbellett
Linux registered user 299762 @ Linux registered machine 213003

 [...]
hmm...
io ti sconsiglierei questa strada dell'inheritance; meglio l'embedding 
in questi casi (IMHO) - cioe': secondo me, usa pure i valarray se vuoi 
per implementare gli elementi della matrice, pero' usa il valarray come 
*membro*:
<CODE>
template <typename T>
class Matrix
{
public:
  ...
  // Ritorna # righe
  int rows() const;
  // Ritorna # colonne
  int cols() const;
  // Accedi ad un elemento della matrice
  T operator( int r, int c ) const;
  T & operator( int r, int c );
  // ...
protected:
  valarray<T> m; // elementi della matrice
  int numRows;   // # righe
  int numCols;   // # colonne
};
</CODE>
In questa maniera, sei anche libero di cambiare l'implementazione della 
classe, es: al posto di un unico valarray usi un valarray di valarray 
(es.: matrice vista come colonna di righe, o riga di colonne), oppure 
usi un'altra classe ancora. Insomma: svincoli la tua classa matrix dalla 
sua implementazione interna.
Se devi definire operator+, lo fai tra matrici, e, in *questa* 
particolare implementazione, sfrutterai + di valarray
nella classe Matrix:
<CODE>
friend Matrix operator+( const Matrix & lhs, const Matrix & rhs );
</CODE>
Quando implementi operator+:
<CODE>
// ...corpo di operaotr+
{
    // ...fai i controlli sulle dimensioni delle matrici in input
    // Suppongo che Matrix abbia un costruttore che prenda le
    // dimensioni della matrice
    Matrix S( righe, colonne ); // crei matrice somma
    // Sommi i valarray :)
    S.m = lhs.m + rhs.m;
    // Ritorna copia di S
    return S;
}
</CODE>
Piccola nota: l'uso di reference count per gestire le matrici potrebbe 
essere importante se le matrici sono di grosse dimensioni, e vuoi 
evitare costoste chiamate a costruttori di copia.
Ciao,
Dan

 Dan
infatti ... se faccio somme ed operazioni su valarray (solo con alcuni 
compilatori).. va veloce uguale (se non addirittura di piu') 
dell'equivalente elemento per elemento...
questa tecnica del ref counting mi interessa parecchio... si puo' 
applicare anche a funzioni generiche che definisco io, non membro?
una persona del ng mi aveva suggerito una specie di breve tutorial in 
rete sul ref counting ma non ci ho capito molto. Sapresti suggerirmi 
qualche referenza (anche cartacea) ?
grazie


 
ehmmm ... dicesi pattern proxy?


 qualche referenza (anche cartacea) ?
:( purtroppo no... Credo pero' che se cerchi con Google troverai un 
sacco di esempi e tutorial!
...a dire il vero, c'era un vecchio articolo di Computer Programming, 
probabilmente tradotto da Dr. Dobb's, in cui si parlava di ref. counting 
proprio applicati ad una classe per gestire i vettori.
...ma e' un articolo di fine anni 90 probabilmente ...(ne ho un vago 
ricordo).
Ciao,
Dan
PS: Vado in vacanza per una decina di giorni, non credo quindi che 
potro' frequentare il NG da li'; se hai altre domande ti rispondero' al 
ritorno.

 ritorno.
ho trovato in rete una libreria di smart pointers
http://ootips.org/yonat/4dev/smart-pointers.html#What
uno dei quali implementa il copy on write con reference counting....
   template <class T>
   class matrix  {
       size_t dim1_,dim2_;
       cow_ptr<valarray <T> > dati_;
}
dove cow_ptr č lo smart pointer che dovrebbe garantirmi l'efficienza...
tuttavia questa efficienza non č proprio raggiunta... forse perchč 
sbaglio a fare qualcosa nella redifinizione dell'operatore +
matrix <T> inline operator+(const matrix <T> &input1,const  matrix <T> 
&input2)
{
   cow_ptr<valarray <T> > risultato(new valarray <T> 
(*(input1.dati_)+*(input2.dati_)));
   return matrix<T>(input1.dim1_,input1.dim2_,risultato); }


 
L'uso di 'cow_ptr' evita la copia del 'valarray' quando viene restituito 
l'oggetto 'matrix', ma non impedisce la creazione dell'oggetto temporaneo 
restituito da *(input1.dati_)+*(input2.dati_).
Potresti fare così:
  valarray <T>*  tmp = new valarray <T> (input1.dati_ );
  *tmp += input2.dati_;
  cow_ptr< valarray <T> > risultato( tmp );
  return matrix<T>(input1.dim1_,input1.dim2_,risultato);
Max

 Max
e comunque non cambia niente ...
ho provato ad aggiungere un cout << "copy"<<endl; nella funzione copy di 
cow_ptr e effettivamente copia solo alla fine delle addizioni...
boh ... sono arrivato alla conclusione che l'unica maniera per 
utilizzare oggetti grandi con prestazioni portabili senza farsi venire 
mal di testa.. è quella di fare le operazioni elemento per elemento e di 
fare dei brutti ma efficaci
void funzione(const &paperino,const &pluto,&pippo)
per simulare la funzione
pippo=funzione(paperino,pluto)
alla fine il buon principio del KISS vince
la mia classetta derivata da valarray offre delle prestazioni elevate e 
molto portabili per quanto riguarda la creazione e l'accesso agli 
elementi in un file di poche righe che mi posso portare ovunque 
appresso.. alla fine a me serve questo... il resto è solo curiosità e 
sfizio di imparare tecniche nuove
mi piacerebbe capire cosa hanno combinato quelli della gnu per la classe 
valarray (ho provato a leggere il codice sorgente... un casino 
incomprensibile) per renderla cosi' performante...
è come se quando scrivo
pippo=pluto+ minnie + topolino
il compilatore automaticamente traduca tutto in
for (size_t i=0;i<pippo.size();i++) pippo[i]=pluto[i]+minnie[i]+topolino[i];
P.S. Gli smart pointer della boost sono migliori di quelli che ho 
scaricato dal sito?


basandomi su un contenitore prelevato dalla  TNT ho finalmente raggiunto 
il mio scopo... scrivere una classe di matrici efficiente e semplice, 
con la quale ottenere prestazioni portabili con i vari compilatori, 
indipendentemente da quanto sia ottimizzata la locale valarray.
I contenitori forniti dalla TNT, sia per interfaccia, sia per la 
struttura di contenitori di contenitori (progettata per ottenere la 
compatibilità di notazione con gli array multidimensionali del C) non mi 
soddisfacevano.
fmatrix e farray sono indicizzate a partire da 1, fmatrix è 
column-oriented ed è pronta da usare con lapack e blas.
per esempio data la funzione sstev_
extern "C"
{
   void sstev_(char &job, int &n, float *d, float *e,  float *z, int 
&ldz, float *work,int &info);
}
che prende i vettori d, z e work e la matrice "e" basta chiamarla
come
sstev_ (job,n, &d[0], &e[0], &z[0], ldz, &work[0],info);
Ho provato ad implementare la metaprogrammazione con template per 
espandere operazioni del tipo A=B+C-D*(E/F) pero' era talmente 
farraginoso da ottenere che dopo essermi scervellato per ottenere 
qualcosa che funzionasse solo con l'addizione ho preferito desistere.
In ogni caso l'operatore [] 0-based permette di scorrere in maniera 
naturale tutti gli elementi del contenitore qualsiasi tipo esso sia con 
la formula
for (size_t i=0;size_t<pippo.size();i++) a[i]=b[i]+c[i]-d[i]*(e[i]/f[i])
con il vantaggio che questa formula è indipendente da come sono fatti a, 
b ,c , d ,e,f ... l'importante è che abbiano lo stesso numero di 
elementi (matrici 10*10, array fortran di 100 elementi, matrici fortran 
5*20 vanno tutte bene)
/*  -*-C++-*-
These classes are provided by Antonio Suriano as his own expansion of TNT.
  I really wanted a simple ref counted copy on write class. I do not use 
template metaprogramming to expand series of operations like A=B+C+D ... 
better to use
manual expansion ... i provide a generic [] to let people easily make 
this by a for (size_t i=0;i<A.size();i++) A[i]=B[i]+C[i]+D[i]
Do not use copy creator entering function because of references... 
create a new object and return this.
Function calling is zero cost.
class i_refvec is provided by Pozo and belongs to TNT project
*/
/*
*
* Template Numerical Toolkit (TNT)
*
* Mathematical and Computational Sciences Division
* National Institute of Technology,
* Gaithersburg, MD USA
*
*
* This software was developed at the National Institute of Standards and
* Technology (NIST) by employees of the Federal Government in the course
* of their official duties. Pursuant to title 17 Section 105 of the
* United States Code, this software is not subject to copyright protection
* and is in the public domain. NIST assumes no responsibility whatsoever for
* its use by other parties, and makes no guarantees, expressed or implied,
* about its quality, reliability, or any other characteristic.
*
*/
/*
Modified and expanded by Antonio Suriano
*/
#include <cmath>
using std::sqrt;
#ifndef TNT_I_REFVEC_H
#define TNT_I_REFVEC_H
#include <cstdlib>
#include <iostream>
#ifdef TNT_BOUNDS_CHECK
#include <cassert>
#endif
#ifndef NULL
#define NULL 0
#endif
namespace tonymatrix
{
/*
	Internal representation of ref-counted array.  The TNT
	arrays all use this building block.
	<p>
	If an array block is created by TNT, then every time
	an assignment is made, the left-hand-side reference
	is decreased by one, and the right-hand-side refernce
	count is increased by one.  If the array block was
	external to TNT, the refernce count is a NULL pointer
	regardless of how many references are made, since the
	memory is not freed by TNT.
 
*/
template <class T>
class i_refvec
{
   private:
     T* data_;
     int *ref_count_;
   public:
			 i_refvec();
	explicit i_refvec(size_t n);
	inline	 i_refvec(T* data);
	inline	 i_refvec(const i_refvec &v);
	inline   T*		 begin();
	inline const T* begin() const;
	inline  T& operator[](size_t i);
	inline const T& operator[](size_t i) const;
	inline  i_refvec<T> & operator=(const i_refvec<T> &V);
		    void copy_(T* p, const T* q, const T* e);
		    void set_(T* p, const T* b, const T* e);
	inline 	int	 ref_count() const;
	inline  int is_null() const;
	inline  void destroy();
			 ~i_refvec();
 
};
template <class T>
void i_refvec<T>::copy_(T* p, const T* q, const T* e)
{
	for (T* t=p; q<e; t++, q++)
		*t= *q;
}
template <class T>
i_refvec<T>::i_refvec() : data_(NULL), ref_count_(NULL) {}
/**
	In case n is 0 or negative, it does NOT call new.
*/
template <class T>
i_refvec<T>::i_refvec(size_t n) : data_(NULL), ref_count_(NULL)
{
	if (n >= 1)
	{
#ifdef TNT_DEBUG
		std::cout  << "new data storage.\n";
#endif
		data_ = new T[n];
		ref_count_ = new int;
		*ref_count_ = 1;
	}
}
template <class T>
inline	 i_refvec<T>::i_refvec(const i_refvec<T> &V): data_(V.data_),
	ref_count_(V.ref_count_)
{
	if (V.ref_count_ != NULL)
	    (*(V.ref_count_))++;
}
template <class T>
i_refvec<T>::i_refvec(T* data) : data_(data), ref_count_(NULL) {}
template <class T>
inline T* i_refvec<T>::begin()
{
	return data_;
}
template <class T>
inline const T& i_refvec<T>::operator[](size_t i) const
{
	return data_[i];
}
template <class T>
inline T& i_refvec<T>::operator[](size_t i)
{
	return data_[i];
}
template <class T>
inline const T* i_refvec<T>::begin() const
{
	return data_;
}
template <class T>
i_refvec<T> & i_refvec<T>::operator=(const i_refvec<T> &V)
{
	if (this == &V)
		return *this;
	if (ref_count_ != NULL)
	{
		(*ref_count_) --;
		if ((*ref_count_) == 0)
			destroy();
	}
	data_ = V.data_;
	ref_count_ = V.ref_count_;
	if (V.ref_count_ != NULL)
	    (*(V.ref_count_))++;
	return *this;
}
template <class T>
void i_refvec<T>::destroy()
{
	if (ref_count_ != NULL)
	{
#ifdef TNT_DEBUG
		std::cout << "destorying data... \n";
#endif
		delete ref_count_;
#ifdef TNT_DEBUG
		std::cout << "deleted ref_count_ ...\n";
#endif
		if (data_ != NULL)
			delete []data_;
#ifdef TNT_DEBUG
		std::cout << "deleted data_[] ...\n";
#endif
		data_ = NULL;
	}
}
/*
* return 1 is vector is empty, 0 otherwise
*
* if is_null() is false and ref_count() is 0, then
*
*/
template<class T>
int i_refvec<T>::is_null() const
{
	return (data == NULL ? 1 : 0);
}
/*
*  returns -1 if data is external,
*  returns 0 if a is NULL array,
*  otherwise returns the positive number of vectors sharing
*  		this data space.
*/
template <class T>
int i_refvec<T>::ref_count() const
{
	if (data_ == NULL)
		return 0;
	else
		return (ref_count_ != NULL ? *ref_count_ : -1) ;
}
template <class T>
i_refvec<T>::~i_refvec()
{
	if (ref_count_ != NULL)
	{
		(*ref_count_)--;
		if (*ref_count_ == 0)
		destroy();
	}
}
#endif
/* TNT_I_REFVEC_H */
template <class T> class matrix;
template <class T> class array;
template <class T> class fmatrix;
template <class T> class farray;
   template <class T> matrix <T>operator+(const matrix <T> &input1,const 
matrix <T> &input2);
   template <class T> matrix <T>operator-(const matrix <T> &input1,const 
matrix <T> &input2);
   template <class T> matrix <T>operator*(const matrix <T> &input1,const 
matrix <T> &input2);
   template <class T> matrix <T>operator/(const matrix <T> &input1,const 
matrix <T> &input2);
   template <class T> fmatrix <T>operator+(const fmatrix <T> 
&input1,const fmatrix <T> &input2);
   template <class T> fmatrix <T>operator-(const fmatrix <T> 
&input1,const fmatrix <T> &input2);
   template <class T> fmatrix <T>operator*(const fmatrix <T> 
&input1,const fmatrix <T> &input2);
   template <class T> fmatrix <T>operator/(const fmatrix <T> 
&input1,const fmatrix <T> &input2);
   template <class T> array <T>operator+(const array <T> &input1,const 
array <T> &input2);
   template <class T> array <T>operator-(const array <T> &input1,const 
array <T> &input2);
   template <class T> array <T>operator*(const array <T> &input1,const 
array <T> &input2);
   template <class T> array <T>operator/(const array <T> &input1,const 
array <T> &input2);
   template <class T> farray <T>operator+(const farray <T> &input1,const 
farray <T> &input2);
   template <class T> farray <T>operator-(const farray <T> &input1,const 
farray <T> &input2);
   template <class T> farray <T>operator*(const farray <T> &input1,const 
farray <T> &input2);
   template <class T> farray <T>operator/(const farray <T> &input1,const 
farray <T> &input2);
template <class T>
class matrix
{
  private:
   i_refvec <T> dati_;
   size_t dim1_,dim2_;
  public:
   explicit matrix (const size_t a,const size_t b) :dim1_(a), dim2_(b), 
dati_(a*b) {}
   matrix (const size_t a,const size_t b,const  i_refvec <T>  input) 
:dim1_(a), dim2_(b),dati_(input) {}
   T& operator() (const size_t i,const  size_t j) {return dati_[i*dim2_+j];}
   const  T& operator() (const size_t i,const size_t j) const {return 
dati_[i*dim2_+j];}
   T& operator[] (const size_t i){return dati_[i];}
   const  T& operator[] (const size_t i) const {return dati_[i];}
   size_t dim1() const {return dim1_;}
   size_t dim2() const {return dim2_;}
   size_t size() const {return dim1_*dim2_;};
   matrix copy()const  {matrix res(dim1_,dim2_); for (size_t 
i=0;i<size();i++) res[i]=dati_[i]; return res;}
   matrix operator+=(const matrix<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]+=input.dati_[i]; return *this;}
   matrix operator-=(const matrix<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]-=input.dati_[i]; return *this;}
   matrix operator*=(const matrix<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]*=input.dati_[i]; return *this;}
   matrix operator/=(const matrix<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]/=input.dati_[i]; return *this;}
   matrix operator+=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]+=input; return *this;}
   matrix operator-=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]-=input; return *this;}
   matrix operator*=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]*=input; return *this;}
   matrix operator/=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]/=input; return *this;}
   friend matrix operator+<>(const matrix <T> &input1,const matrix <T> 
&input2);
   friend matrix operator-<>(const matrix <T> &input1,const matrix <T> 
&input2);
   friend matrix operator*<>(const matrix <T> &input1,const matrix <T> 
&input2);
   friend matrix operator/<>(const matrix <T> &input1,const matrix <T> 
&input2);
   matrix (const matrix<T> &input): dim1_(input.dim1_), 
dim2_(input.dim2_), dati_(input.dati_) {}
   T dot_product(const matrix <T> &input) const {T res=0; for (size_t 
i=0;i<size();i++) res+=dati_[i]*input.dati_[i]; return res;}
   T norm() const {return sqrt(dot_product(*this));}
   void normalize(){operator/=(norm());}
   void set(const  T input){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=input;}
   void apply( T func(T)){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=func(dati_[i]);}
   void apply( T func()){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=func();}
};
  template <class T>
  matrix <T> inline operator+(const matrix <T> &input1,const matrix <T> 
&input2)
  {
    matrix <T> risultato(input2.dim1_,input2.dim2_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]+input2[i];
    return risultato;
  }
  template <class T>
  matrix <T> inline operator-(const matrix <T> &input1,const matrix <T> 
&input2)
  {
    matrix <T> risultato(input2.dim1_,input2.dim2_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]-input2[i];
    return risultato;
  }
  template <class T>
  matrix <T> inline operator*(const matrix <T> &input1,const matrix <T> 
&input2)
  {
    matrix <T> risultato(input2.dim1_,input2.dim2_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]*input2[i];
    return risultato;
  }
  template <class T>
  matrix <T> inline operator/(const matrix <T> &input1,const matrix <T> 
&input2)
  {
    matrix <T> risultato(input2.dim1_,input2.dim2_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]/input2[i];
    return risultato;
  }
template <class T>
class fmatrix
{
  private:
   i_refvec <T> dati_;
   size_t dim1_,dim2_;
  public:
     explicit fmatrix (const size_t a,const size_t b) :dim1_(a), 
dim2_(b), dati_(a*b) {}
     fmatrix (const size_t a,const size_t b,const  i_refvec <T>  input) 
:dim1_(a), dim2_(b),dati_(input) {}
     T& operator() (const size_t i,const  size_t j) {return 
dati_[j*dim1_+j-1-dim1_];}
     const  T& operator() (const size_t i,const size_t j) const {return 
dati_[j*dim1_+j-1-dim1_];}
     T& operator[] (const size_t i){return dati_[i];}
     const  T& operator[] (const size_t i) const {return dati_[i];}
     size_t dim1() const {return dim1_;}
     size_t dim2() const {return dim2_;}
     size_t size() const {return dim1_*dim2_;};
     fmatrix copy()const  {fmatrix res(dim1_,dim2_); for (size_t 
i=0;i<size();i++) res[i]=dati_[i]; return res;}
     fmatrix operator+=(const fmatrix<T> &input){size_t size_=size(); 
for (size_t i=0;i<size_;i++) dati_[i]+=input.dati_[i]; return *this;}
     fmatrix operator-=(const fmatrix<T> &input){size_t size_=size(); 
for (size_t i=0;i<size_;i++) dati_[i]-=input.dati_[i]; return *this;}
     fmatrix operator*=(const fmatrix<T> &input){size_t size_=size(); 
for (size_t i=0;i<size_;i++) dati_[i]*=input.dati_[i]; return *this;}
     fmatrix operator/=(const fmatrix<T> &input){size_t size_=size(); 
for (size_t i=0;i<size_;i++) dati_[i]/=input.dati_[i]; return *this;}
     fmatrix operator+=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]+=input; return *this;}
     fmatrix operator-=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]-=input; return *this;}
     fmatrix operator*=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]*=input; return *this;}
     fmatrix operator/=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]/=input; return *this;}
     friend fmatrix operator+<>(const fmatrix <T> &input1,const fmatrix 
<T> &input2);
     friend fmatrix operator-<>(const fmatrix <T> &input1,const fmatrix 
<T> &input2);
     friend fmatrix operator*<>(const fmatrix <T> &input1,const fmatrix 
<T> &input2);
     friend fmatrix operator/<>(const fmatrix <T> &input1,const fmatrix 
<T> &input2);
     fmatrix (const fmatrix<T> &input): dim1_(input.dim1_), 
dim2_(input.dim2_), dati_(input.dati_) {}
   T dot_product(const fmatrix <T> &input) const {T res=0; for (size_t 
i=0;i<size();i++) res+=dati_[i]*input.dati_[i]; return res;}
   T norm() const {return sqrt(dot_product(*this));}
   void normalize(){operator/=(norm());}
   void set(const  T &input){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=input;}
   void apply( T func(T)){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=func(dati_[i]);}
   void apply( T func()){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=func();}
};
  template <class T>
  fmatrix <T> inline operator+(const fmatrix <T> &input1,const fmatrix 
<T> &input2)
  {
    fmatrix <T> risultato(input2.dim1_,input2.dim2_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]+input2[i];
    return risultato;
  }
  template <class T>
  fmatrix <T> inline operator-(const fmatrix <T> &input1,const fmatrix 
<T> &input2)
  {
    fmatrix <T> risultato(input2.dim1_,input2.dim2_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]-input2[i];
    return risultato;
  }
  template <class T>
  fmatrix <T> inline operator*(const fmatrix <T> &input1,const fmatrix 
<T> &input2)
  {
    fmatrix <T> risultato(input2.dim1_,input2.dim2_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]*input2[i];
    return risultato;
  }
  template <class T>
  fmatrix <T> inline operator/(const fmatrix <T> &input1,const fmatrix 
<T> &input2)
  {
    fmatrix <T> risultato(input2.dim1_,input2.dim2_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]/input2[i];
    return risultato;
  }
template <class T>
class farray
{
  private:
   i_refvec <T> dati_;
   size_t dim1_;
  public:
     explicit farray (const size_t a) :dim1_(a),  dati_(a) {}
     farray (const size_t a,const  i_refvec <T>  input) :dim1_(a), 
dati_(input) {}
     T& operator() (const size_t i) {return dati_[i-1];}
     const  T& operator() (const size_t i) const {return dati_[i-1];}
     T& operator[] (const size_t i){return dati_[i];}
     const  T& operator[] (const size_t i) const {return dati_[i];}
     size_t dim1() const {return dim1_;}
     size_t size() const {return dim1_;};
     farray copy()const  {farray res(dim1_); for (size_t 
i=0;i<size();i++) res[i]=dati_[i]; return res;}
     farray operator+=(const farray<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]+=input.dati_[i]; return *this;}
     farray operator-=(const farray<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]-=input.dati_[i]; return *this;}
     farray operator*=(const farray<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]*=input.dati_[i]; return *this;}
     farray operator/=(const farray<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]/=input.dati_[i]; return *this;}
     farray operator+=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]+=input; return *this;}
     farray operator-=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]-=input; return *this;}
     farray operator*=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]*=input; return *this;}
     farray operator/=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]/=input; return *this;}
     friend farray operator+<>(const farray <T> &input1,const farray <T> 
&input2);
     friend farray operator-<>(const farray <T> &input1,const farray <T> 
&input2);
     friend farray operator*<>(const farray <T> &input1,const farray <T> 
&input2);
     friend farray operator/<>(const farray <T> &input1,const farray <T> 
&input2);
     farray (const farray<T> &input): dim1_(input.dim1_), 
dati_(input.dati_) {}
   T dot_product(const farray <T> &input) const {T res=0; for (size_t 
i=0;i<size();i++) res+=dati_[i]*input.dati_[i]; return res;}
   T norm() const {return sqrt(dot_product(*this));}
   void normalize(){operator/=(norm());}
   void set( const T &input){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=input;}
   void apply( T func(T)){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=func(dati_[i]);}
   void apply( T func()){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=func();}
};
  template <class T>
  farray <T> inline operator+(const farray <T> &input1,const farray <T> 
&input2)
  {
    farray <T> risultato(input2.dim1_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]+input2[i];
    return risultato;
  }
  template <class T>
  farray <T> inline operator-(const farray <T> &input1,const farray <T> 
&input2)
  {
    farray <T> risultato(input2.dim1_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]-input2[i];
    return risultato;
  }
  template <class T>
  farray <T> inline operator*(const farray <T> &input1,const farray <T> 
&input2)
  {
    farray <T> risultato(input2.dim1_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]*input2[i];
    return risultato;
  }
  template <class T>
  farray <T> inline operator/(const farray <T> &input1,const farray <T> 
&input2)
  {
    farray <T> risultato(input2.dim1_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]/input2[i];
    return risultato;
  }
template <class T>
class array
{
  private:
   i_refvec <T> dati_;
   size_t dim1_;
  public:
     explicit array (const size_t a) :dim1_(a),  dati_(a) {}
     array (const size_t a,const  i_refvec <T>  input) :dim1_(a), 
dati_(input) {}
     T& operator() (const size_t i) {return dati_[i-1];}
     const  T& operator() (const size_t i) const {return dati_[i];}
     T& operator[] (const size_t i){return dati_[i];}
     const  T& operator[] (const size_t i) const {return dati_[i];}
     size_t dim1() const {return dim1_;}
     size_t size() const {return dim1_;};
     array copy()const  {array res(dim1_); for (size_t i=0;i<size();i++) 
res[i]=dati_[i]; return res;}
     array operator+=(const array<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]+=input.dati_[i]; return *this;}
     array operator-=(const array<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]-=input.dati_[i]; return *this;}
     array operator*=(const array<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]*=input.dati_[i]; return *this;}
     array operator/=(const array<T> &input){size_t size_=size(); for 
(size_t i=0;i<size_;i++) dati_[i]/=input.dati_[i]; return *this;}
     array operator+=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]+=input; return *this;}
     array operator-=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]-=input; return *this;}
     array operator*=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]*=input; return *this;}
     array operator/=(const T &input){size_t size_=size(); for (size_t 
i=0;i<size_;i++) dati_[i]/=input; return *this;}
     friend array operator+<>(const array <T> &input1,const array <T> 
&input2);
     friend array operator-<>(const array <T> &input1,const array <T> 
&input2);
     friend array operator*<>(const array <T> &input1,const array <T> 
&input2);
     friend array operator/<>(const array <T> &input1,const array <T> 
&input2);
     array (const array<T> &input): dim1_(input.dim1_), 
dati_(input.dati_) {}
   T dot_product(const array <T> &input) const {T res=0; for (size_t 
i=0;i<size();i++) res+=dati_[i]*input.dati_[i]; return res;}
   T norm() const {return sqrt(dot_product(*this));}
   void normalize(){operator/=(norm());}
   void set( const T &input){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=input;}
   void apply( T func(T)){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=func(dati_[i]);}
   void apply( T func()){size_t a=size();for (size_t i=0;i<a;i++) 
dati_[i]=func();}
};
  template <class T>
  array <T> inline operator+(const array <T> &input1,const array <T> 
&input2)
  {
    array <T> risultato(input2.dim1_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]+input2[i];
    return risultato;
  }
  template <class T>
  array <T> inline operator-(const array <T> &input1,const array <T> 
&input2)
  {
    array <T> risultato(input2.dim1_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]-input2[i];
    return risultato;
  }
  template <class T>
  array <T> inline operator*(const array <T> &input1,const array <T> 
&input2)
  {
    array <T> risultato(input2.dim1_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]*input2[i];
    return risultato;
  }
  template <class T>
  array <T> inline operator/(const array <T> &input1,const array <T> 
&input2)
  {
    array <T> risultato(input2.dim1_);
    size_t size=input1.size();
    for (size_t i=0;i<size;i++) risultato[i]=input1[i]/input2[i];
    return risultato;
  }
} /* namespace tonymatrix */
