/[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 1486 by artak, Thu Apr 10 00:10:50 2008 UTC revision 1488 by artak, Fri Apr 11 00:22:31 2008 UTC
# Line 583  def GMRES(b, Aprod, Msolve, bilinearform Line 583  def GMRES(b, Aprod, Msolve, bilinearform
583     m=iter_restart     m=iter_restart
584     iter=0     iter=0
585     while True:     while True:
586        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
587        x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)        x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
588        iter+=iter_restart        iter+=iter_restart
589        if stopped: break        if stopped: break
# Line 597  def GMRESm(b, Aprod, Msolve, bilinearfor Line 597  def GMRESm(b, Aprod, Msolve, bilinearfor
597     norm_b=math.sqrt(r_dot_r)     norm_b=math.sqrt(r_dot_r)
598    
599     if x==None:     if x==None:
600        x=0*b        x=0
601     else:     else:
602        r=Msolve(b-Aprod(x))        r=Msolve(b-Aprod(x))
603        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
# Line 610  def GMRESm(b, Aprod, Msolve, bilinearfor Line 610  def GMRESm(b, Aprod, Msolve, bilinearfor
610     v=[]     v=[]
611    
612     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
613      
614     v.append(r/rho)     v.append(r/rho)
615     g[0]=rho     g[0]=rho
616    
# Line 683  def GMRESm(b, Aprod, Msolve, bilinearfor Line 684  def GMRESm(b, Aprod, Msolve, bilinearfor
684     else : xhat=v[0]     else : xhat=v[0]
685            
686     x += xhat     x += xhat
687     if iter!=iter_restart-1:     if iter<iter_restart-1:
688        stopped=True        stopped=True
689     else:     else:
690        stopped=False        stopped=False
# Line 1051  class HomogeneousSaddlePointProblem(obje Line 1052  class HomogeneousSaddlePointProblem(obje
1052        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1053                return self.__reduction*self.getTolerance()                return self.__reduction*self.getTolerance()
1054    
1055        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG'):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=10):
1056                """                """
1057                solves the saddle point problem using initial guesses v and p.                solves the saddle point problem using initial guesses v and p.
1058    
# Line 1080  class HomogeneousSaddlePointProblem(obje Line 1081  class HomogeneousSaddlePointProblem(obje
1081                self.iter=0                self.iter=0
1082            if solver=='GMRES':                  if solver=='GMRES':      
1083                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1084                  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.,iter_restart=iter_restart)
1085                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1086                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1087                  #       u=v+(u-v)                  #       u=v+(u-v)
# Line 1093  class HomogeneousSaddlePointProblem(obje Line 1094  class HomogeneousSaddlePointProblem(obje
1094                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1095                  #       u=v+(u-v)                  #       u=v+(u-v)
1096          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1097                else:                
1098                  if solver=='PCG':
1099                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1100                  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)
1101              u=r[0]                u=r[0]  

Legend:
Removed from v.1486  
changed lines
  Added in v.1488

  ViewVC Help
Powered by ViewVC 1.1.26