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

revision 1514 by artak, Wed Apr 16 00:15:44 2008 UTC revision 1517 by artak, Fri Apr 18 02:36:37 2008 UTC
# Line 758  def MINRES(b, Aprod, Msolve, bilinearfor Line 758  def MINRES(b, Aprod, Msolve, bilinearfor
758      #---------------------------------------------------------------------      #---------------------------------------------------------------------
759      # Main iteration loop.      # Main iteration loop.
760      # --------------------------------------------------------------------      # --------------------------------------------------------------------
761      while not stoppingcriterium(rnorm,Anorm*ynorm):    #  checks ||r|| < (||A|| ||x||) * TOL      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
762
763      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
764          iter    = iter  +  1          iter    = iter  +  1
# Line 888  def TFQMR(b, Aprod, Msolve, bilinearform Line 888  def TFQMR(b, Aprod, Msolve, bilinearform
888  #  #
889  #  while ( iter < kmax-1 ):  #  while ( iter < kmax-1 ):
890
891    while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b):    while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
892      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
893
894      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
# Line 944  def TFQMR(b, Aprod, Msolve, bilinearform Line 944  def TFQMR(b, Aprod, Msolve, bilinearform
944
945      iter = iter + 1      iter = iter + 1
946
print iter
947    return x    return x
948
949
1139        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
1140            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
1141
1142        def __stoppingcriterium_GMRES(self,norm_r,norm_b):        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES'):
1143            return self.stoppingcriterium_GMRES(norm_r,norm_b)            return self.stoppingcriterium2(norm_r,norm_b,solver)

def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
return self.stoppingcriterium_MINRES(norm_r,norm_Ax)

1144
1145        def setTolerance(self,tolerance=1.e-8):        def setTolerance(self,tolerance=1.e-8):
1146                self.__tol=tolerance                self.__tol=tolerance
1151        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1152                return self.__reduction*self.getTolerance()                return self.__reduction*self.getTolerance()
1153
1154        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=10):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1155                """                """
1156                solves the saddle point problem using initial guesses v and p.                solves the saddle point problem using initial guesses v and p.
1157
1180                self.iter=0                self.iter=0
1181            if solver=='GMRES':                  if solver=='GMRES':
1182                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1183                  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)                  p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1184                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1185                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1186                  #       u=v+(u-v)                  #       u=v+(u-v)
1187          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1188
1189            if solver=='TFQMR':                  if solver=='TFQMR':
1190                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1191                  p=TFQMR(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)                  p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1192                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1193                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1194                  #       u=v+(u-v)                  #       u=v+(u-v)
1196
1197            if solver=='MINRES':                  if solver=='MINRES':
1198                  if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1199                  p=MINRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)                  p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1200                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1201                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1202                  #       u=v+(u-v)                  #       u=v+(u-v)
1214        def __Msolve(self,r):        def __Msolve(self,r):
1215            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1216
1217        def __Msolve_GMRES(self,r):        def __Msolve2(self,r):
1218            return self.solve_prec(r)            return self.solve_prec(r)
1219
1220