/************************************************************************

  Copyright
  Alessandro MIRONE
  mirone@esrf.fr

  Copyright 2002  by European Synchrotron Radiation Facility, Grenoble, 
                  France

                               ----------
 
                           All Rights Reserved
 
                               ----------

Permission to use, copy, modify, and distribute this software and its
documentation for any purpose and without fee is hereby granted,
provided that the above copyright notice appear in all copies and that
both that copyright notice and this permission notice appear in
supporting documentation, and that the names of European Synchrotron
Radiation Facility or ESRF or BLISS not be used in advertising or 
publicity pertaining to distribution of the software without specific, 
written prior permission.

EUROPEAN SYNCHROTRON RADIATION FACILITY DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL EUROPEAN SYNCHROTRON
RADIATION FACILITY OR ESRF BE LIABLE FOR ANY SPECIAL, INDIRECT OR 
CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, 
DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.

**************************************************************************/
/***************************************************************************
                          Sparsa.h  -  description
                             -------------------
    begin                : Tue Feb 1 2000
    copyright            : (C) 2000 by Alessandro MIRONE
    email                : mirone@lure.u-psud.fr
 ***************************************************************************/
#ifndef _SPARSAH_
#define _SPARSAH_
class MatriceSparsa {
public:
	MatriceSparsa() {threads=0;};
	virtual void preparathreads(int t) {};
	virtual ~MatriceSparsa();
	virtual void Moltiplica(double *ris, double *vect  );
	virtual void solveCGSym( double *x, double *b, double tol, int MAXIT, double *initial=0);	
	virtual void trasforma(double fattore, double addendo);
	virtual void copiasu( MatriceSparsa * &m);
  int dim;
  int dim2;
  int threads;
};

class DoubleProvider {
public:
  DoubleProvider() {};
  DoubleProvider(int nsteps, int dim){};
  
  void puliscimemoria(int nsteps, int dim){};
  virtual void settamemoria(int nsteps, int dim){};
  virtual double * operator [] (int i ){return NULL;};
  virtual ~DoubleProvider() {};
};

class DoubleProviderSimple: public DoubleProvider {
public:

  DoubleProviderSimple();
  DoubleProviderSimple(int nsteps, int dim);
  void puliscimemoria();

  void settamemoria(int nsteps, int dim);
  double * operator [] (int i );
  
  ~DoubleProviderSimple();

 private:
  double **q;
  int nsteps;
  int dim;
};

class DoubleCache: public DoubleProvider {
public:
  
  DoubleCache();
  DoubleCache(int nsteps, int dim);
  void puliscimemoria();

  void settamemoria(int nsteps, int dim);
  double * operator [] (int i );


  ~DoubleCache();
private:
  double *q, * qref;


  int accesscounter; // da incrementare ogni volta che si fa []
  int icached[10]; // da inizializzare a -1
  int axcached[10];
  int  magazzino;
  int dim;
  int nsteps;
  char nametemplate[200];
};




class LanczosSteps;


class LanczosSteps{
public:
	LanczosSteps();
	inline void setMatrice(MatriceSparsa *m, int  ns, int tipo=0 ) { matrice=m; inizializza(); assicuramemoria(matrice->dim, ns, tipo); };
	inline void setInversione(double t, int m) { inversione=1; tol=t; maxit=m;};
	inline void setDiretto      () { inversione=0;};
	void GramSchmidt(double* z, int n ) ;
	void GramSchmidt( int n ) ;
	void passeggia(int k, int m,  double * iniziale, int continued=0   );
	void passeggiaDebug(int k, int m,  double * iniziale, int continued=0   );
	void inizializza();
	~LanczosSteps(	);

    friend class Lanczos;
    void setNoOldDiags() { olddiags=0;};
    void diagonalizza(int k, int m);
    int convergeds(int m, double tole=2.0e-16);
    void ricipolla(int k, int m);
    int lineaestnulla(int k, double tole);
public :
	void liberamemoria();
	void assicuramemoria(int d, int ns, int tipo=0);
	MatriceSparsa * matrice;
	int dim;

	int nsteps;
	
	// double **q;
	DoubleProvider *q;
	double *alpha;
	double *oldalpha;
	int olddiags;
	double * beta;
	
	double *p;
	double **omega;

	int inversione;
	double tol;
	int        maxit;
	
	uMatrix<double> A;
	uMatrix<double> eval;
	uMatrix<double>  evect;

	int compact;
};

class Lanczos{
public:
	Lanczos();
	inline void setMatrice(MatriceSparsa *m, int inversione=0, double t=1.0e-16 , int max=40000)
	 {
	 		matrice=m;
	 		if(inversione==1) {
	 			steps.setInversione(t,max);
	 		} else  {
	 			steps.setDiretto();
	 		}	 		
	 };
  ~Lanczos(	);
	int cerca(int nd, double shift=0.0, double tol=2.0e-16);
	int cercaDebug(int nd, double shift=0.0, double tol=2.0e-16);
	double * geteigenvector(int n) { return (* (steps.q))[n];}
	double    &geteigenvalue(int n)  { return steps.eval(n,n); };
	LanczosSteps  steps;
private :
	// void liberamemoria();
	// void assicuramemoria(int d, int ns);
	
	MatriceSparsa * matrice;
	MatriceSparsa * shiftata;
	int nd;
	int nsteps;
	
	int inversione;
	double shift;
	double tol;
	int        maxit;
};



class SparsaS;

class Array ;

class Sparsa3A : public  MatriceSparsa {

 public:
  Sparsa3A();
  Sparsa3A(int N, double *c, int *i, int *j);
   ~Sparsa3A();
  void inizializza(char *file, double fact=1);
  void somma(char *file , double f=1 ) ;
  void sommaFast(char *file , double f=1 ) ;

  void inizializza( Sparsa3A &a ) ;
  void inizializza( Sparsa3A &rea,Sparsa3A &ima  ) ;
  void inizializza( Sparsa3A &a , double fact) ;
  void inizializza_bin(char *file  , double fact=1.0) ;

  void somma    ( Sparsa3A &a, double f=1 ) ;

  void pulisci();
  void Moltiplica(double *ris, double *vect  );
  void MoltiplicaMinus(double *ris, double *vect  );
  void MoltiplicaDiag     (double *ris, double *vect  );
  void MoltiplicaDiagMinus(double *ris, double *vect  );



  void Moltiplica(Array *ris, Array *vect  );
  void MoltiplicaMinus(Array *ris, Array *vect  );

  void MoltiplicaDiag(Array * ris, Array*vect  );
  void MoltiplicaDiagMinus(Array * ris, Array*vect  );

  void Moltiplica_i(double *ris, double *vect  );

  void aggiungiElemento    (double a, int i, int j);
  void aggiungiElementoFast(double a, int i, int j);

  void riordina();
  void trim();
  void solveCGSym( double *x, double *b, double tol, int MAXIT, double *initial=0);
  void solveCGSym( double *x,double *y,double *b, double *d, double tol, int MAXIT, double *initial=0);
  void solveCGSym2menoImag( double *x, double *b, double imag, double tol, int MAXIT, double *initial=0);
  void trasforma(double fattore, double addendo) ;	
  void copiasu( MatriceSparsa * & m);
  void scrivi(char *file) ;
  void scrivi() ;
  void transpose() { int *dum; dum=col; col=row; row=dum; };
  void raddoppia(double imag) ;
  void preparathreads(int t);
  void Moltiplica(double *ris, double *vect , int start, int end );

  int * intervallithreads;
  friend class SparsaS;

  void   gohersch   ();
  double goherschMin();
  double goherschMax();

  double getelement(int i,int j);

	void solveBCG( double *x,double *b, double tol, int itmax,int & iter,double &err);

	void ottieniuM(uMatrix<double> &a);

// private:

  void reqmrpbCG( double *r, double *b);
  void aggiungimemoria();
	double dAd(double *d, double *dum );
	double dA2menoImag(double *d, double *dum , double imag, double * preco);
	std::complex<double>  dAd(double *d0, double *di0, double *d, double *dum , double *dumi);

  int n;
  int nsize;

  double *coeff  ;
  int    *col   ;
  int    *row   ;

  double *Gmin, *Gmax;
  int nG;
};

void dumptofile(int N,  double *ArrayFLOAT ,  char * nome ) ;
void dumptofile(double *iniziale, int dim, char * nome ) ;

#include<stdexcept>


struct Array {
  Array( int size ) { 
    firstdimension=0;
    ownsdata=1;
    init_array( size ); 
  }

  Array( int fd, int sd ) { 
    firstdimension=fd;
    ownsdata=1;
    init_array( fd*sd ); 
  }
  Array( int size , double * data) { 
    firstdimension=0;
    ownsdata=0;
    this->data=data;
    this->size=size;
  }

  ~Array() { 
    if( ownsdata) 
      free_array(); 
    else
      std::cout << " cancellata una finestra su array \n";
  }


  Array * get_array(int index) {
    if( index >=firstdimension || index<0) {
      throw std::out_of_range("index invalido in Array * get_array(int index)");
    }
    return new Array(  size/firstdimension, data+ index*(size/firstdimension)  );
  }
  
  
  double get_value( int index ) const
  {
    printf("get_value(%d)\n",index);
    return data[index];
  }
  
  void set_value( int index, double value )
  {
    printf("set_value(%d,%e)\n",index,value);
    data[index]=value;
  }
  
#undef RAND_MAX
#define RAND_MAX        2147483647
  
  void set_all_random(  double value )
  {
    
    int index;
    for(index= 0; index<size; index++) {
      data[index]=value*random()*1.0/RAND_MAX;
    }
    
  }
  void set_to_zero(  )
  {
    int index;
    for(index= 0; index<size; index++) {
      data[index]=0.0;
    }      
  }
  void set_to_one(  )
  {
    int index;
    for(index= 0; index<size; index++) {
      data[index]=1.0;
    }      
  }
  
  static void copy_to_a_from_b(Array &a, Array&b) {
    if( a.size != b.size) {
      throw std::out_of_range("frames are not aligned");   
    }
    memcpy(a.data, b.data, a.size*sizeof(double) );
  }
  
  
  void copy_from_b( Array&b) {
    if( this->size != b.size) {
      throw std::out_of_range("frames are not aligned");   
    }
    memcpy(this->data, b.data, this->size*sizeof(double) );
  }
  
  
  void add_from_vect( Array&vect) {
    if( size != vect.size) {
      throw std::out_of_range("frames are not aligned");   
    }
    for(int i=0;i<size;++i)
      data[i]+=vect.data[i];
  }
  
  void add_from_vect_with_fact( Array&vect, double fact) {
    if( size != vect.size) {
      throw std::out_of_range("frames are not aligned");   
    }
    for(int i=0;i<size;++i)
      data[i]+=vect.data[i]*fact;
  }
  void mult_by_fact(  double fact){
    for(int i=0;i<size;++i)
      data[i]*=fact;
  }
  
  
  void  set_indices_inv(int N, int *indices, double *values) ;
  void  set_indices( int  N, int *indices, double *values) ;
  
  
  Array* copy() const
  {
    Array *b=new Array(size);
    b->firstdimension=firstdimension;
    for(int i=0;i<size;++i)
      b->data[i]=data[i];
    return b;
  }
  
  Array& operator+=( const Array& a ) {
    for(int i=0;i<len() && i<a.len();++i) data[i]+=a.data[i];
    return *this; 
  }
  int len() const { return size; }
  
  double *dataAddress() {return data;}
  
  void dumptofile(char * nome) {
    FILE * f=fopen(nome,"w");
    fwrite(dataAddress(), sizeof(double), len(), f);
    fclose(f);
    
  }
  void  dividebyarray( Array &b);
  static double scalare(Array &a, Array &b);
  static double sqrtscalare(Array &a, Array &b);
  static void mat_mult(Array & b, double * ArrayFLOAT2dNOCHECK ,Array & a);



  void normalizza(double norm);
  void normalizzaauto();
  
protected:

  int firstdimension;
  int ownsdata;
  int size;
  double* data;
  
  void init_array( int n )
  {
    ownsdata=1;
    size=n;
    data=(double*)malloc(sizeof(double)*n);
    memset(data, 0, sizeof(double)*n);
    printf("array %p, size=%d, data=%p\n", this, size, data);
  }
  void free_array()
  {
    printf("free_array: %p\n", this );
    free(data);
  }
  
};

typedef Array ListArray;


















class Sparsa2A : public  MatriceSparsa {

 public:
  Sparsa2A(int N, double *c, int *i, int *j, int dim, int dim2);
  ~Sparsa2A();

  void inizializza( Sparsa2A &a ) ;

  void pulisci();
  void Moltiplica(double *ris, double *vect  );

  void solveCGSym( double *x, double *b, double tol, int MAXIT, double *initial=0){}; 

    void trasforma(double fattore, double addendo){} ;	

    void copiasu( MatriceSparsa * & m){};

  friend class SparsaS;

  void   gohersch   ();
  double goherschMin();
  double goherschMax();

  double getelement(int i,int j);
  
//  private:
  
  void reqmrpbCG( double *r, double *b);

  double dAd(double *d, double *dum );
  double dA2menoImag(double *d, double *dum , double imag, double * preco);

  int n;
  int nsize;

  int dim, dim2;

  double *coeff  ;
  int    *row   ;
  int *ends;

  double *Gmin, *Gmax;
  int nG;
};











class SparsaS : public  MatriceSparsa  {

 public:
  
  SparsaS();
 SparsaS( double *x, double *b, double tol, int MAXIT, double *initial=0);
 virtual ~SparsaS();

  void inizializza(Sparsa3A &s);
  void inizializza(SparsaS   &s);

  void Moltiplica(double *ris, double *vect  );

  void copiasu( MatriceSparsa * & m);

  friend class Sparsa3A;
 private:

  void resize(int np, int nhorsp );
  int n;
  int nhors;
  double *coeff  ;
  int    *ncol   ;
};



class Sparsa_Resolver {
public:
  void resize(int n);
  void azzera();
  void sommaH(int  i,int j, double add);
  void diagonalizzaH ( double Dene, uMatrix<double>& eVect, uMatrix<double>& eVal , int maxNevals=-1, double tole=2.0e-15,
                        int minDim=5, double distFromPole=0.001 );
  double prodottoscalare(  uMatrix<double>& a , uMatrix<double>&    b    );
  inline void stampa(){ 
    matricione.scrivi();
   }
  void salvasufile(char *nome);
  void addizionadafile(char *nome, double f);
 inline void somma(Sparsa_Resolver & a, double f=1) { matricione.somma(a.matricione,f); }
 void inizializza(Sparsa3A & m);
private:
  Sparsa3A   matricione;
};





void diagonalizza4py(int k, int m,double *alpha,  double *beta, double *pevect, double *peval);
double sqrtscalare(double * a, double *b, int n);
double scalare(double * , double *, int);
double scalare( int , double * , double *);
void normalizza(double * , double , int);
void normalizzaauto(double * ,  int);
void somma      (double * ,    double * ,    double , int);

#define NO_SPARSAMPI
#ifdef SPARSAMPI
#define MAXNSPARSA3AP 1000
#define MAXNVECTORP   10000
class inizializzaMPI {
 public:
  inizializzaMPI();
  ~inizializzaMPI();

};
class Sparsa3AP: public  MatriceSparsa {
 public:

  Sparsa3AP();
  ~Sparsa3AP();

  void caricaArrays(int Nels, int * from, int * to, double * coeffs,int dim_from, int dim_to);
  void riceviArrays();


  int item_id;
  int proc_id;
  int n_procs;

  static Sparsa3AP** items;
 private:
  void crea();

  int dim_from;
  int dim_to;
  int *Nelspp;
  int **topp;
  int **frompp;
  double ** coeffspp;

  int *dims2;
  int ** loc2glob;

  int **others_loc2glob;
  int *others_dims2;

};



class VectorP {
 public:
  VectorP(int Nels);
  ~VectorP();
  void caricaArray(int Nels,  double * coeffs);
  void riceviArrays();
  void setToZero();
  void add( double factor, VectorP * v   );
  void scale( double factor  );

  friend double scalarP(VectorP *a, VectorP *b);



  int item_id;
  int proc_id;
  int n_procs;
  

  static VectorP** items;
 private:
  void crea(int Nels);
  int dim;
  double *coeffs;
};

double scalarP(VectorP *a, VectorP *b);





#endif

#endif

