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

  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 Array ;

class Sparsa3A  {

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

  void inizializza( Sparsa3A &a ) ;
  void inizializza( Sparsa3A &a , double fact) ;
  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 trasforma(double fattore, double addendo) ;	
  void transpose() { int *dum; dum=col; col=row; row=dum; };
  void Moltiplica(double *ris, double *vect , int start, int end );

  void aggiungimemoria();

  int n;
  int nsize;
  int dim;
  int dim2;


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

};



#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;


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

