/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Diff of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1466 by artak, Wed Apr 2 03:28:25 2008 UTC revision 1467 by gross, Wed Apr 2 08:10:37 2008 UTC
# Line 49  import util Line 49  import util
49  import math  import math
50    
51  ##### Added by Artak  ##### Added by Artak
52  from Numeric import zeros,Int,Float32,Float64  # from Numeric import zeros,Int,Float64
53  ###################################  ###################################
54    
55    
# Line 551  type like argument C{x}. Line 551  type like argument C{x}.
551    
552  #Apply a sequence of k Givens rotations, used within gmres codes  #Apply a sequence of k Givens rotations, used within gmres codes
553  # vrot=givapp(c, s, vin, k)  # vrot=givapp(c, s, vin, k)
554  def givapp(c,s,vin,k):  def givapp(c,s,vin):
555      vrot=vin      vrot=vin # warning: vin is altered!!!!
556      for i in range(k+1):      if isinstance(c,float):
557          w1=c[i]*vrot[i]-s[i]*vrot[i+1]          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
558      w2=s[i]*vrot[i]+c[i]*vrot[i+1]      else:
559          vrot[i:i+2]=w1,w2          for i in range(len(c)):
560                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
561            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
562                vrot[i:i+2]=w1,w2
563      return vrot      return vrot
564    
565  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
   
    from numarray import dot  
     
    v0=b  
566     iter=0     iter=0
567       r=Msolve(b)
568       r_dot_r = bilinearform(r, r)
569       if r_dot_r<0: raise NegativeNorm,"negative norm."
570       norm_b=math.sqrt(r_dot_r)
571    
572     if x==None:     if x==None:
573        x=0*b        x=0*b
574     else:     else:
575        b += (-1.)*Aprod(x)        r=Msolve(b-Aprod(x))
576     r=b        r_dot_r = bilinearform(r, r)
577          if r_dot_r<0: raise NegativeNorm,"negative norm."
    rhat=Msolve(r)  
     
    rhat_dot_r = bilinearform(rhat, rhat)  
    norm_r=math.sqrt(rhat_dot_r)  
   
    if rhat_dot_r<0: raise NegativeNorm,"negative norm."  
578        
579     h=zeros((iter_max,iter_max),Float32)     h=numarray.zeros((iter_max,iter_max),numarray.Float64)
580     c=zeros(iter_max,Float32)     c=numarray.zeros(iter_max,numarray.Float64)
581     s=zeros(iter_max,Float32)     s=numarray.zeros(iter_max,numarray.Float64)
582     g=zeros(iter_max,Float32)     g=numarray.zeros(iter_max,numarray.Float64)
583     v=[]     v=[]
584    
585     v.append(rhat/norm_r)     rho=math.sqrt(r_dot_r)
586     rho=norm_r     v.append(r/rho)
587     g[0]=rho     g[0]=rho
588    
589     while not stoppingcriterium(rho,rho,norm_r):     while not stoppingcriterium(rho,norm_b):
590    
591      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
592    
593            
594      vhat=Aprod(v[iter])      p=Msolve(Aprod(v[iter]))
     p=Msolve(vhat)    
595    
596      v.append(p)      v.append(p)
597    
# Line 627  def GMRES(b, Aprod, Msolve, bilinearform Line 624  def GMRES(b, Aprod, Msolve, bilinearform
624      if iter > 0 :      if iter > 0 :
625          hhat=[]          hhat=[]
626          for i in range(iter+1) : hhat.append(h[i][iter])          for i in range(iter+1) : hhat.append(h[i][iter])
627          hhat=givapp(c[0:iter],s[0:iter],hhat,iter-1);          hhat=givapp(c[0:iter],s[0:iter],hhat);
628              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i][iter]=hhat[i]
629    
630      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
# Line 636  def GMRES(b, Aprod, Msolve, bilinearform Line 633  def GMRES(b, Aprod, Msolve, bilinearform
633          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1][iter]/mu
634          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
635          h[iter+1][iter]=0.0          h[iter+1][iter]=0.0
636          g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2],0)          g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
637    
638  # Update the residual norm  # Update the residual norm
639          rho=abs(g[iter+1])          rho=abs(g[iter+1])
# Line 646  def GMRES(b, Aprod, Msolve, bilinearform Line 643  def GMRES(b, Aprod, Msolve, bilinearform
643  # It's time to compute x and leave.          # It's time to compute x and leave.        
644    
645     if iter > 0 :     if iter > 0 :
646       y=zeros(iter,Float32)         y=numarray.zeros(iter,numarray.Float64)    
647       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1][iter-1]
648       if iter > 1 :         if iter > 1 :  
649          i=iter-2            i=iter-2  
650          while i>=0 :          while i>=0 :
651            y[i] = ( g[i] - dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
652            i=i-1            i=i-1
653       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
654       for i in range(iter-1):       for i in range(iter-1):
# Line 854  class HomogeneousSaddlePointProblem(obje Line 851  class HomogeneousSaddlePointProblem(obje
851        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
852            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
853    
854        def __stoppingcriterium_GMRES(self,norm_r,rho,r):        def __stoppingcriterium_GMRES(self,norm_r,norm_b):
855            return self.stoppingcriterium_GMRES(rho,r)            return self.stoppingcriterium_GMRES(norm_r,norm_b)
856    
857        def setTolerance(self,tolerance=1.e-8):        def setTolerance(self,tolerance=1.e-8):
858                self.__tol=tolerance                self.__tol=tolerance
# Line 892  class HomogeneousSaddlePointProblem(obje Line 889  class HomogeneousSaddlePointProblem(obje
889                self.iter=0                self.iter=0
890            if solver=='GMRES':                  if solver=='GMRES':      
891                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
892                  p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1)                  p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
893                    # solve Au=f-B^*p
894                    #       A(u-v)=f-B^*p-Av
895                    #       u=v+(u-v)
896          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
897            
898                else:                else:
899                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
900                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p*1)                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
901              u=r[0]                u=r[0]  
                 u=v+self.solve_A(v,p)     # Lutz to check !!!!!  
902    
903            return u,p            return u,p
904    
# Line 914  class HomogeneousSaddlePointProblem(obje Line 913  class HomogeneousSaddlePointProblem(obje
913            # return BA^-1B*p            # return BA^-1B*p
914            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
915            v=self.solve_A(self.__z,p)            v=self.solve_A(self.__z,p)
916            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, -self.B(v))
917    
918        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
919            # return BA^-1B*p            # return BA^-1B*p
920            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
921        v=self.solve_A(self.__z,p)        v=self.solve_A(self.__z,p)
922            return self.B(v)            return -self.B(v)
923    
924  class SaddlePointProblem(object):  class SaddlePointProblem(object):
925     """     """
Legend:
Removed from v.1466  
changed lines
  Added in v.1467
  ViewVC Help
Powered by ViewVC 1.1.26