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

revision 1469 by gross, Thu Apr 3 05:16:56 2008 UTC revision 1475 by artak, Mon Apr 7 05:13:21 2008 UTC
# Line 579  def givapp(c,s,vin): Line 579  def givapp(c,s,vin):
579              vrot[i:i+2]=w1,w2              vrot[i:i+2]=w1,w2
580      return vrot      return vrot
581
582  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
583       m=iter_restart
584       iter=0
585       while True:
586          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
587          x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
588          iter+=iter_restart
589          if stopped: break
590       return x
591
592    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
593     iter=0     iter=0
594     r=Msolve(b)     r=Msolve(b)
595     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
# Line 593  def GMRES(b, Aprod, Msolve, bilinearform Line 603  def GMRES(b, Aprod, Msolve, bilinearform
603        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
604        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm,"negative norm."
605
606     h=numarray.zeros((iter_max,iter_max),numarray.Float64)     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
607     c=numarray.zeros(iter_max,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
608     s=numarray.zeros(iter_max,numarray.Float64)     s=numarray.zeros(iter_restart,numarray.Float64)
609     g=numarray.zeros(iter_max,numarray.Float64)     g=numarray.zeros(iter_restart,numarray.Float64)
610     v=[]     v=[]
611
612     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
613     v.append(r/rho)     v.append(r/rho)
614     g[0]=rho     g[0]=rho
615
616     while not stoppingcriterium(rho,norm_b):     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
617
618      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
619
# Line 673  def GMRES(b, Aprod, Msolve, bilinearform Line 683  def GMRES(b, Aprod, Msolve, bilinearform
683     else : xhat=v[0]     else : xhat=v[0]
684
685     x += xhat     x += xhat
686       if iter!=iter_restart-1:
687          stopped=True
688       else:
689          stopped=False
690
691     return x     return x,stopped
692
693  #############################################  #############################################
694
# Line 919  class HomogeneousSaddlePointProblem(obje Line 933  class HomogeneousSaddlePointProblem(obje
933                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
934                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
935              u=r[0]                u=r[0]
936                print "div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
937
938            return u,p            return u,p
939

Legend:
 Removed from v.1469 changed lines Added in v.1475