import math
import profile
import time
import Hilbertxx
import Sparsa
import Numeric

######################################
## interfaccia per sparse
# __class__
# static    .getClass4Vect()
# object    .gohersch()
# object    .goherschMax()
# object    .goherschMin()

######################################
## interfaccia per vettori
# __class__ === class4vect
# class Vector(dim)
# object.set_value(n,v)
# object.set_all_random(v)

def REAL(a):
    if( type(a)==type(1.0+1.0j) ):
        return a.real
    else:
        return a

class Callable:
    def __init__(self, anycallable):
        self.__call__=anycallable






class TrialVectorC:

    def __init__(self, *dim):
        if( dim!=(0,) ):
            self.vr=Sparsa.Array(*dim)
            self.vi=Sparsa.Array(*dim)
        else:
            pass
        
    def __getitem__(self, i):
        res=TrialVectorC(0)
        res.vr=self.vr[i]
        res.vi=self.vi[i]
        return res

    def mat_mult(self, evect , q):
        Sparsa.Array.mat_mult(self.vr, evect , q.vr)
        Sparsa.Array.mat_mult(self.vi, evect , q.vi)

    def dividebyarray(self,prec):
        self.vr.dividebyarray(prec)
        self.vi.dividebyarray(prec)


    def len(self):
	return self.vr.len()

    def copy(self):
	duma=self.vr.copy()
	dumb=self.vi.copy()
        res=TrialVectorC(0)
        res.vr=duma
        res.vi=dumb
        return res 

    def copy_to_a_from_b(self,b):
        Sparsa.Array.copy_to_a_from_b(self.vr, b.vr  )
        Sparsa.Array.copy_to_a_from_b(self.vi, b.vi  )

    def set_value(self, n, val):
        if( type(val)==type(1+1j)):
            self.vr.set_value(n, val.real)
            self.vi.set_value(n, val.imag)
        else:
            self.vr.set_value(n, val)
            self.vi.set_value(n, 0.0)

    def set_to_zero(self):
        self.vr.set_to_zero()
        self.vi.set_to_zero()
        
    def set_to_one(self):
        self.vr.set_to_one()
        self.vi.set_to_zero()
        
    def set_all_random(self, v):
        self.vr.set_all_random(v)
        self.vi.set_all_random(v)
        
    
    def scalare(self,b         ):
        resR =  Sparsa.Array.scalare(self.vr, b.vr)+ Sparsa.Array.scalare(self.vi, b.vi)
        resI =  Sparsa.Array.scalare(self.vr, b.vi)- Sparsa.Array.scalare(self.vi, b.vr)
        return resR + resI * 1.0j

    def sqrtscalare(self,b):
        assert( self is b)
        return math.sqrt( self.scalare(b).real )

    def normalizzaauto(self):
        norma = self.sqrtscalare(self)
        self.vr.normalizza(norma)
        self.vi.normalizza(norma)
        

    def normalizza(self,norma):
        self.vr.normalizza(norma)
        self.vi.normalizza(norma)

    def rescale(self,fact):
      if(fact==0.0):
        self.vr.set_to_zero()
        self.vi.set_to_zero()
      else:
        if type(fact)==type(1.0+1.0j):
          norma=1.0/fact.real
        else:
          norma=1.0/fact
        self.vr.normalizza(norma)
        self.vi.normalizza(norma)
        
    def add_from_vect(self, b):
            self.vr.add_from_vect(b.vr   )
            self.vi.add_from_vect(b.vi     )
        
     
        
    def add_from_vect_with_fact(self, b, fact):
        if( type(fact)==type(1.0+1.0j)):
            self.vr.add_from_vect_with_fact(b.vr  ,    fact.real   )
            self.vi.add_from_vect_with_fact(b.vi  ,    fact.real   )
            self.vr.add_from_vect_with_fact(b.vi  ,    -fact.imag   )
            self.vi.add_from_vect_with_fact(b.vr  ,    fact.imag   )
        else:
            self.vr.add_from_vect_with_fact(b.vr  ,    fact   )
            self.vi.add_from_vect_with_fact(b.vi  ,    fact   )
        
class TrialMatrixC:
    def __init__(self, mR, mI):
        self.mR=mR
        self.mI=mI
        self.dim= self.mR.dim
        self.shift=0.0

    def Moltiplica(self,res,v):
        self.mR.Moltiplica(res.vr,v.vr)
        self.mR.Moltiplica(res.vi,v.vi)
        
        if( self.mI is not None):
          self.mI.Moltiplica(res.vi,v.vr)
          self.mI.MoltiplicaMinus(res.vr,v.vi)

        if( self.shift !=0.0):
          res.add_from_vect_with_fact(v,self.shift)

    def MoltiplicaMinus(self,res,v):
        self.mR.MoltiplicaMinus(res.vr,v.vr)
        self.mR.MoltiplicaMinus(res.vi,v.vi)
        
        if( self.mI is not None):
          self.mI.MoltiplicaMinus(res.vi,v.vr)
          self.mI.Moltiplica(res.vr,v.vi)

        if( self.shift !=0.0):
          res.add_from_vect_with_fact(v,-self.shift)

    def MoltiplicaDiag(self,res,v):
        self.mR.MoltiplicaDiag(res.vr,v.vr)
        self.mR.MoltiplicaDiag(res.vi,v.vi)
        
        if( self.mI is not None):
          self.mI.MoltiplicaDiag(res.vi,v.vr)
          self.mI.MoltiplicaDiagMinus(res.vr,v.vi)

        if( self.shift !=0.0):
          res.add_from_vect_with_fact(v,self.shift)

    def MoltiplicaDiagMinus(self,res,v):
        self.mR.MoltiplicaDiagMinus(res.vr,v.vr)
        self.mR.MoltiplicaDiagMinus(res.vi,v.vi)
        
        if( self.mI is not None):
          self.mI.MoltiplicaDiagMinus(res.vi,v.vr)
          self.mI.MoltiplicaDiag     (res.vr,v.vi)

        if( self.shift !=0.0):
          res.add_from_vect_with_fact(v,self.shift)
         
        
    def gohersch(self):
        self.mR.gohersch()
        print " dim R est " , self.mR.dim
        print " dim I est " , self.mI.dim
        if( self.mI is not None):
          self.mI.gohersch()
        
    def goherschMax(self):
      if( self.mI is not None):
        return self.mR.goherschMax()+self.mI.goherschMax()
      else:
        return self.mR.goherschMax()

    def trasforma(self, *args):
        self.mR.trasforma(*args)

                   
    def getClass4Vect():
        return TrialVectorC
    getClass4Vect=Callable(getClass4Vect)



class TrialMatrix(Sparsa.Sparsa3A):
    def getClass4Vect():
        return Sparsa.Array
    getClass4Vect=Callable(getClass4Vect)

def stampa(z):
     for i in range(5):
         print z.vr.get_value(i)


def biCGsolve(A, At,  b,  tol=1.0e-10, nmax = 1000, prec=None, restart=None):
    dim=A.dim
    print dim
    
    class4sparse = A.__class__
    class4vect   = class4sparse.getClass4Vect()
    
    r=class4vect(dim)
    rbar=class4vect(dim)
    
    w=class4vect(dim)
    wbar=class4vect(dim)

    x = class4vect(dim)



    if restart is None:
      x.set_to_zero()
    else:
      x.copy_to_a_from_b(restart)

    
    p=class4vect(dim)
    pbar=class4vect(dim)
    
    dum=class4vect(dim)
    z    = class4vect(dim)
    zbar = class4vect(dim)
    
    # r =  b - dot(self, x)
    
    r.copy_to_a_from_b(b)
    A.MoltiplicaMinus(r,x)

    rbar.copy_to_a_from_b(r) # manca il coniugato che fovrebbe essere preso da scalar

    


    z   .copy_to_a_from_b(r)
    if not (prec is  None):
        z.dividebyarray(prec)



    
    k = 0
    err=10
    while abs(err) > tol and k < nmax:




        # zz = rr/kvect
        zbar.copy_to_a_from_b(rbar)
        if not (prec is None):        
            zbar.dividebyarray(prec)


       
        rho = class4vect.scalare(rbar,z).real

        print rho
        if( abs(rho)<1.0e-15):
          return x

        beta=0
        if(k==0):
            p.copy_to_a_from_b(z)

            pbar.copy_to_a_from_b(zbar)

        else:
            beta=rho/rhoold
            
            p.rescale(beta)
            p.add_from_vect(z)
            
            pbar.rescale(beta)
            pbar.add_from_vect(zbar)
            
        rhoold=rho

        z.set_to_zero()
        A.Moltiplica(z,p)
         
        alpha=(rho/(class4vect.scalare(pbar,z).real))



        
        zbar.set_to_zero()
        At.Moltiplica(zbar,pbar)


        x.add_from_vect_with_fact(p,alpha)
        r   .add_from_vect_with_fact(z   ,-alpha)
        rbar.add_from_vect_with_fact(zbar,-alpha)




        class4vect.scalare(pbar,z).real


        z.copy_to_a_from_b(r)


        if prec is not None:
            z.dividebyarray(prec)

        dum.set_to_zero()
        A.Moltiplica(dum,x)
        dum.add_from_vect_with_fact(b,-1)
        err=class4vect.scalare(dum,dum)
        
        print k,"  ", err
        k = k+1
        
    return x








class TrialMatrix(Sparsa.Sparsa3A):
    def getClass4Vect():
        return Sparsa.Array
    getClass4Vect=Callable(getClass4Vect)

    
class Lanczos:
    dump_count=0;
    countdumpab=0
    def __init__(self, sparse):

        self.matrice=sparse

        self.class4sparse = sparse.__class__
        self.class4vect   = self.class4sparse.getClass4Vect()


        self.nsteps = 0
        self.dim=self.matrice.dim

        self.q = None

        self.alpha = None
        self.beta = None
        self.omega=None
                
        self.tol = 1e-15
        self.maxIt = 50

        self.old = False
        



    def diagoCustom(self, minDim=5):

        self.matrice.gohersch()
        shift = - self.matrice.goherschMax()
        self.cerca(minDim, shift)

        for ne in range(len(self.eval)):
            self.eval[ne] =  self.eval[ne] - shift

    def passeggia(self, k, m, start, gram=0):

        if k<0 or m>self.nsteps:
            print "k = ",k," m = ",m," nsteps = ",self.nsteps
            print "Something wrong in passe k<0 or m>nsteps"
            sys.exit(1)

        sn   = math.sqrt(float(self.dim))
        eu   = 1.1e-16
        eusn = eu*sn
        eq   = math.sqrt(eu)


        if k==0:
            
            
            self.class4vect.copy_to_a_from_b(self.q[0],start)
            
            self.q[0].normalizzaauto( )

            # self.q[0].dumptofile("q_0")



        p= self.class4vect(self.dim)

        for i in range(k,m):
            p.set_to_zero()
            # self.q[i].dumptofile("qbef_"+str(self.dump_count) )
           
            self.matrice.Moltiplica(p,self.q[i])

            self.alpha[i] = REAL(self.class4vect.scalare(p , self.q[i]))
            
            # p.dumptofile("p"+str(self.dump_count) )
            self.dump_count+=1

            
            if i==k:
                p.add_from_vect_with_fact( self.q[k] ,   -self.alpha[k]     ) 
                for l in range(k):
                    p.add_from_vect_with_fact( self.q[l] ,  -self.beta[l]     ) 

            else:
                # self.q[i].dumptofile("q_i"+str(self.dump_count) )
                # self.q[i-1].dumptofile("q_im1"+str(self.dump_count) )
                p.add_from_vect_with_fact(self.q[i]  ,    -self.alpha[i]    ) 
                p.add_from_vect_with_fact(self.q[i-1]  ,  -self.beta[i-1]   ) 


            # p.dumptofile("pp"+str(self.dump_count) )
            self.dump_count+=1

            
            self.beta[i]=self.class4vect.sqrtscalare(p,p) 

            self.omega[i,i]=1.

            max0 = 0.0

            if self.beta[i] != 0:  #and  not scipy.isnan(self.beta[i]):
                for j in range(i+1):
                    self.omega[i+1,j] = eusn



                    if j<k:
                        add = 2 * eusn + abs(self.alpha[j]-self.alpha[i]) \
                              * abs(self.omega[i,j])
                        if i!=k:
                            add += self.beta[j]*abs(self.omega[i,k])
                        if i>0 and j!=i-1:
                            add += self.beta[i-1]*abs(self.omega[i-1,j])
                        self.omega[i+1,j] += add / self.beta[i]

                       
                    elif j==k :                        
                        add = 2 * eusn + abs(self.alpha[j]-self.alpha[i]) \
                              * abs(self.omega[i,j])

                        for w in range(k):
                            add += self.beta[w]*abs(self.omega[i,w])

                        if i!=k+1:
                            add += self.beta[k]*abs(self.omega[i,k+1])
                        if i>0 and i!=k+1:
                            add += self.beta[i-1]*abs(self.omega[i-1,k])

                        self.omega[i+1,j] += add / self.beta[i]



                    elif j<i :

                        add = 2 * eusn + abs(self.alpha[j]-self.alpha[i]) \
                              * abs(self.omega[i,j])

                        if i!=j+1:
                            add += self.beta[j]*abs(self.omega[i,j+1])

                        if i>0 and j>0:
                            add += self.beta[j-1]*abs(self.omega[i-1,j-1])

                        if i>0 and i!=j+1:
                            add += self.beta[i-1]*abs(self.omega[i-1,j])

                        self.omega[i+1,j] += add / self.beta[i]


                            
                    else:

                        add = eusn

                        if i>0:
                            add += self.beta[i-1]*abs(self.omega[i,i-1] )

                        self.omega[i+1,j] += add / self.beta[i]


                    self.omega[j,i+1] = self.omega[i+1,j]

                    max0 += self.omega[i+1,j]**2
            if self.beta[i]==0 or max0>eu or gram :
##            if self.beta[i]==0 or max0>0. :

                if i>0:
                    print " GRAMO  self.beta[i]==0 or max0>eu and i>0", i,"  ", self.dump_count
                    self.GramSchmidt(self.q[i],i)

                    self.q[i].normalizzaauto()

                    p.set_to_zero()
                    
                    self.matrice.Moltiplica(p, self.q[i])

                    self.alpha[i] =  REAL(self.class4vect.scalare(p,self.q[i]))

                print " GRAMO bis ", i
                self.GramSchmidt(p,i+1)
                
                 
                self.beta[i] = self.class4vect.sqrtscalare(p,p)
                p .normalizzaauto( )

                if i>0:
                    condition = eu * \
                                math.sqrt(self.dim * (self.alpha[i]**2+\
                                                      self.beta[i-1]**2))
                else:
                    condition = eu * math.sqrt(self.dim * (self.alpha[i]**2))


                if self.beta[i]< condition:

                    print "starting with a vector perpendicular"

                    self.beta[i]=0.
                    
                    # p.set_to_one()
                    # p.set_value(i,i)
                    p.set_all_random(1.0)

                    
                    self.GramSchmidt(p,i+1)

                    p .normalizzaauto()

                for l in range(i) :
                    self.omega[i,l]=self.omega[l,i]=eusn
                for l in range(i+1):
                    self.omega[i+1,l]=self.omega[l,i+1]=eusn
                    
            else:
                # print " NORMALIZZO PRIMA DI PASSARE ", self.dump_count
                # p.dumptofile("p_bef_norm"+str(self.dump_count) )
                p.normalizzaauto( )
                # p.dumptofile("normprelude"+str(self.dump_count))

            self.class4vect.copy_to_a_from_b(self.q[i+1],p)

          

    def GramSchmidt(self, vect , n):
        for h in range(4):
            for i in range(n):
                s=self.class4vect.scalare( self.q[i], vect)
                vect.add_from_vect_with_fact(self.q[i],-s)


    def allocaMemory(self):
        self.alpha = Numeric.zeros(self.nsteps,Numeric.Float64)
        self.beta  = Numeric.zeros(self.nsteps,Numeric.Float64)
        self.omega = Numeric.zeros((self.nsteps+1,self.nsteps+1),Numeric.Float64)
        self.evect = Numeric.zeros((self.nsteps,self.nsteps),Numeric.Float64)
        self.eval  = Numeric.zeros((self.nsteps),Numeric.Float64)
        self.oldalpha = Numeric.zeros((self.nsteps),Numeric.Float64)
        self.q=self.class4vect(self.nsteps+1,self.dim)


      
    def cerca(self, nd, shift):
        
        self.shift=shift

        self.matrice.trasforma(1.0, shift)

        
        m = min(4*nd, self.dim)

        self.nsteps = m

        self.allocaMemory()


##         self.q=[]
##         for i in range(self.nsteps+1):
##             self.q.append(self.class4vect(self.dim) )



        vect_init = self.class4vect(self.dim)
        # vect_init.set_value(0,1.0)
        
        vect_init.set_all_random(1.0)

        k=0
        nc=0
        self.passeggia(k,m,vect_init)
        
        while nc<nd :

            # Sparsa.dumptofile(self.alpha,  "a"+str(self.countdumpab))  
            # Sparsa.dumptofile(self.beta[:-1] ,  "b"+str(self.countdumpab))  

            print " DIAGONALIZZAZIONE "
            self.diago(k,m)


            # Sparsa.dumptofile(Numeric.resize(self.evect,[m*m]) ,  "vect_"+str(self.countdumpab))  

        
            nc = self.converged(m)

            if k and not Numeric.sometrue(abs(self.beta[:k])>self.tol) :
                break

            if (nc+2*nd) >= m:
                k=m-1
            else:
                k=nc+2*nd

            # self.q.dumptofile("qtotbef"+str(self.countdumpab))
            print "KKKKKKKKK  " ,  k 
            self.ricipolla(k,m)
            # self.q[0].dumptofile("q0afte"+str(self.countdumpab))
            self.countdumpab+=1
            
            print " k,m , dim", k, m, self.dim
            # return m # sentinell
        
            self.passeggia(k,m,self.q[m])

        if m==self.dim:
            return m
        else:
            return k


    def converged(self, m):

        for j in range(1,m):
            for i in range(m-j):
                if abs(self.eval[i]) < abs(self.eval[i+1]):
                    self.eval[i],self.eval[i+1]=self.eval[i+1],self.eval[i]
                    self.evect[i:i+2]=self.evect[i+1],self.evect[i]


        # print "eval", self.eval[0:m]
        res = 0
        if self.old:
            while res<m:
                print self.oldalpha[res]
                if (abs(self.eval[res]-self.oldalpha[res])\
                    /abs(self.oldalpha[res])>self.tol):
                    break
                res+=1
        else:
            self.old=True

        self.oldalpha[0:m]=self.eval[0:m]

        return res

    def ricipolla(self, k,m):

        print "k,m",k,m
        self.alpha[0:k]= self.eval[0:k]


        self.beta[0:k] = self.beta[m-1]*self.evect[0:k,m-1]
        print "beta beta",self.beta,self.beta[m-1],self.evect[0:k,m-1]

        # E = self.evect[0:k,0:m].copy()

        a=self.class4vect(k,  self.dim )

        print " LANCIO mat MULT "

        self.class4vect.mat_mult(a, self.evect[0:k,0:m] , self.q)
## void  Array::mat_mult(Array & b, double * mat,Array & a) {
##   int m=a.firstdimension-1;
##   int k=b.firstdimension;
##   int dim = a.size/a.firstdimension;

##   // uMatrix<double> E(m, k,  mat);
## //   uMatrix<double> A(dim, m, a.data );
## //   uMatrix<double> B(dim, k, b.data );
  
##   ::dumptofile(mat, m*k, "bloccovectpy");
   
##   char transa='N', transb='N';
##   double alpha = double(1.0);
##   double beta = double(0.0);
  
##   dgemm_( &transa, &transb, &dim, &k, &m, &alpha, a.data, &dim, 
##            mat, &m, &beta, b.data, &dim );
## }


        

        for i in range(k):
            self.class4vect.copy_to_a_from_b( self.q[i] , a[i]     )

        a=None

        self.class4vect.copy_to_a_from_b(self.q[k]  ,  self.q[m] ) 

 
        o = self.omega[0:m,0:m].copy() 


        o = Numeric.dot(o,Numeric.transpose(self.evect))


        for i in range(k):
            self.omega[i,k]=self.omega[k,i]=o[i,k]

        
        o = Numeric.dot(self.evect,o)

        
        self.omega[0:k,0:k]=o[0:k,0:k]



         
    def diago(self, k, m):

        eval=Numeric.zeros([m,m],Numeric.Float64)

        print " ALPHA  BETA "
        print self.alpha
        print self.beta
        
        Sparsa.diagonalizza4py(k,self.alpha,  self.beta, self.evect, eval)
        for i in range(m):
            self.eval[i]=eval[i,i]







def solveEigenSystem( S_base , nsearchedeigen): 
    lnczs=Lanczos( S_base )
    lnczs.diagoCustom(minDim=nsearchedeigen)
    return lnczs.eval, lnczs.q




def   faispettri( evals, evects   , S_ecci , maxene, Dipos, tospan, getvects=0, gram=0):
  class4vect= S_ecci.__class__.getClass4Vect()
  print " dim " , S_ecci.dim
  dumA = class4vect(S_ecci.dim)
  data4spect={}
  neigtc=0
  dene = evals - evals[0]
  
  while( neigtc< len( evals) and  dene[neigtc] < maxene):
    neigtc+=1

  data4spect["neig"] = neigtc
  data4spect["ndipo"] = len(Dipos)
  data4spect["enes"] = {}
  data4spect["alphas"] = {}
  data4spect["norms"] = {}
  data4spect["storevect"] = {}



  lnczs=Lanczos( S_ecci )
  lnczs.nsteps = tospan
  lnczs.allocaMemory()

  for ival in range(neigtc):
    data4spect["enes"][ival]=evals[ival]
    data4spect["alphas"][ival]={}
    data4spect["norms"][ival]={}
    data4spect["storevect"][ival]={}

    for idip in range(len(Dipos )):

      dumA.set_to_zero()
      print evects[ival] 
      norm_to_remember = class4vect.scalare(evects[ival] ,  evects[ival] )

      print "NORMA " , norm_to_remember

      if(Dipos[idip] is not None):
        Dipos[idip].Moltiplica(dumA, evects[ival]  )
      else:
        dumA.copy_to_a_from_b(evects[ival]  )
      
      norm_to_remember = class4vect.scalare( dumA,  dumA )

      print "NORMA " , norm_to_remember
      data4spect["norms"][ival][idip] = norm_to_remember
      if( abs(norm_to_remember) > 1.0e-20):
        lnczs .passeggia(0, lnczs.nsteps, dumA, gram);

        data4spect["alphas"][ival][idip]=Numeric.array([ (lnczs.alpha[j], lnczs.beta[j]) for j in range(lnczs.nsteps)    ])
        if getvects:
          data4spect["storevect"][ival][idip]=lnczs.q.copy()

      else:
        data4spect["alphas"][ival][idip]=[ (0.0, 0.0) for j in range(lnczs.nsteps)    ]
        

  return data4spect
      



def spectra(alphas, energies):
    AB=Numeric.array(alphas)
    A=AB[:,0]
    B=AB[:,1]
    res=ContFracSpectra( A,B, energies )
    return res
  
def ContFracSpectra( A,B, W ):
    n=len(A)
    res=0*W
    for i in range(n-1, -1, -1):
        if( i!=(n-1) ):
            res= -B[i]*B[i]/res
        res=res +A[i]-W
    return 1./res
     
def getinvfirstcolumns( alphas, ein, gammain):

    print alphas
    print ein
    print gammain

    N=len(alphas)
    detes=Numeric.zeros(len(alphas)+1,"D")
    detes[N]=1
    detes[N-1]= (alphas[N-1][0]-ein- 1.0j *gammain)
    for i in range(N-2,-1,-1):
      detes[i] = detes[i+1] * (alphas[i][0]-ein- 1.0j* gammain)- detes[i+1]*alphas[i][1]*alphas[i][1]
      print detes[i]
   
    fattori=Numeric.zeros(len(alphas),"D")
    prod=1
    
    for j in range(N):
      fattori[j]= prod*detes[j+1] /detes[0]
      prod=-prod*alphas[j][1]*alphas[j][1]


    return fattori

def CGsolve(H,eps2,  b, tol=1.0e-10, nmax = 1000, restart=None):
        """
        Solve self*x = b and return x using the conjugate gradient method
        """
        class4sparse = H.__class__
        class4vect   = class4sparse.getClass4Vect()

        print " in CG " 
        dim = H.dim
        r= class4vect(dim)
        p = class4vect(dim)
        diff = class4vect(dim)
        dum=class4vect(dim)
        z=class4vect(dim)
        x0=class4vect(dim)

        x0.set_to_zero()
        if restart is not None:
          x0.copy_to_a_from_b(restart )

        # kvec = self.diag() # preconditionner
        x = x0 # initial guess


        #  residup 
        dum.set_to_zero()
        H.Moltiplica(dum,x)
        r.set_to_zero()
        H.Moltiplica(r,dum)
        r.add_from_vect_with_fact(x, eps2)
        r.normalizza(-1.0)
        r.add_from_vect_with_fact(b, 1.0)


        p.set_to_zero()
        beta = 0.0;
                        
        rho = class4vect.scalare(r, r).real


        err= rho
                        

                        
        k = 0
        while abs(err) > tol and k < nmax:
          
          #  p = r + beta*p;

          dum.set_to_zero()
          dum.add_from_vect_with_fact(r,1.0)
          dum.add_from_vect_with_fact(p,beta)
          p.set_to_zero()
          p.add_from_vect_with_fact(dum,1.0)
          
          
          # z = dot(self, p);
          dum.set_to_zero()
          H.Moltiplica(dum,p)
          z.set_to_zero()
          H.Moltiplica(z,dum)
          z.add_from_vect_with_fact(p, eps2)



          # alpha = rho/vector.dot(p, z);
          alpha = rho/(class4vect.scalare(p, z).real)
                                
          # r = r - alpha*z;
          r.add_from_vect_with_fact(z,-alpha)

                                
                                
          rhoold = rho;
          #rho = vector.dot(r, w);
          
          rho = class4vect.scalare(r, r).real


                                
          # x = x + alpha*p;
          x.add_from_vect_with_fact(p,alpha)
                                
          beta = rho/rhoold;
                                
          # err = vector.norm(dot(self, x) - b);



          dum.set_to_zero()
          H.Moltiplica(dum,x)
          diff.set_to_zero()
          H.Moltiplica(diff,dum)
          diff.add_from_vect_with_fact(x, eps2)
          diff.normalizza(-1.0)
          diff.add_from_vect_with_fact(b, 1.0)
          

          err= class4vect.scalare(diff,diff).real
          
          
          print err
          print type(err)
          print k,"  ", err
          k = k+1
        return x




