/[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 1488 by artak, Fri Apr 11 00:22:31 2008 UTC revision 1489 by artak, Mon Apr 14 04:29:30 2008 UTC
# Line 474  class IterationHistory(object): Line 474  class IterationHistory(object):
474    
475         """         """
476         self.history.append(norm_r)         self.history.append(norm_r)
477         if self.verbose: print "iter: %s:  inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])         if self.verbose: print "iter: #s:  inner(rhat,r) = #e"#(len(self.history)-1, self.history[-1])
478         return self.history[-1]<=self.tolerance * self.history[0]         return self.history[-1]<=self.tolerance * self.history[0]
479    
480     def stoppingcriterium2(self,norm_r,norm_b):     def stoppingcriterium2(self,norm_r,norm_b):
# Line 491  class IterationHistory(object): Line 491  class IterationHistory(object):
491    
492         """         """
493         self.history.append(norm_r)         self.history.append(norm_r)
494         if self.verbose: print "iter: %s:  norm(r) = %e"%(len(self.history)-1, self.history[-1])         if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])
495         return self.history[-1]<=self.tolerance * norm_b         return self.history[-1]<=self.tolerance * norm_b
496    
497  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
# Line 844  def MINRES(b, Aprod, Msolve, bilinearfor Line 844  def MINRES(b, Aprod, Msolve, bilinearfor
844    
845      return x      return x
846            
847    
848    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
849    
850    # TFQMR solver for linear systems
851    #
852    #
853    # initialization
854    #
855      errtol = math.sqrt(bilinearform(b,b))
856      norm_b=errtol
857      kmax  = iter_max
858      error = []
859    
860      if math.sqrt(bilinearform(x,x)) != 0.0:
861        r = b - Aprod(x)
862      else:
863        r = b
864    
865      r=Msolve(r)
866    
867      u1=0
868      u2=0
869      y1=0
870      y2=0
871    
872      w = r
873      y1 = r
874      iter = 0
875      d = 0
876      
877      v = Msolve(Aprod(y1))
878      u1 = v
879      
880      theta = 0.0;
881      eta = 0.0;
882      tau = math.sqrt(bilinearform(r,r))
883      error = [ error, tau ]
884      rho = tau * tau
885      m=1
886    #
887    #  TFQMR iteration
888    #
889    #  while ( iter < kmax-1 ):
890      
891      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b):
892        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
893    
894        sigma = bilinearform(r,v)
895    
896        if ( sigma == 0.0 ):
897          raise 'TFQMR breakdown, sigma=0'
898        
899    
900        alpha = rho / sigma
901    
902        for j in range(2):
903    #
904    #   Compute y2 and u2 only if you have to
905    #
906          if ( j == 1 ):
907            y2 = y1 - alpha * v
908            u2 = Msolve(Aprod(y2))
909    
910          m = 2 * (iter+1) - 2 + (j+1)
911          if j==0:
912             w = w - alpha * u1
913             d = y1 + ( theta * theta * eta / alpha ) * d
914          if j==1:
915             w = w - alpha * u2
916             d = y2 + ( theta * theta * eta / alpha ) * d
917    
918          theta = math.sqrt(bilinearform(w,w))/ tau
919          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
920          tau = tau * theta * c
921          eta = c * c * alpha
922          x = x + eta * d
923    #
924    #  Try to terminate the iteration at each pass through the loop
925    #
926         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
927         #   error = [ error, tau ]
928         #   total_iters = iter
929         #   break
930          
931    
932        if ( rho == 0.0 ):
933          raise 'TFQMR breakdown, rho=0'
934        
935    
936        rhon = bilinearform(r,w)
937        beta = rhon / rho;
938        rho = rhon;
939        y1 = w + beta * y2;
940        u1 = Msolve(Aprod(y1))
941        v = u1 + beta * ( u2 + beta * v )
942        error = [ error, tau ]
943        total_iters = iter
944        
945        iter = iter + 1
946    
947      print iter
948      return x
949    
950    
951  #############################################  #############################################
952    
953  class ArithmeticTuple(object):  class ArithmeticTuple(object):
# Line 1087  class HomogeneousSaddlePointProblem(obje Line 1191  class HomogeneousSaddlePointProblem(obje
1191                  #       u=v+(u-v)                  #       u=v+(u-v)
1192          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1193    
1194              if solver=='TFQMR':      
1195                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1196                    p=TFQMR(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
1197                    # solve Au=f-B^*p
1198                    #       A(u-v)=f-B^*p-Av
1199                    #       u=v+(u-v)
1200            u=v+self.solve_A(v,p)
1201    
1202            if solver=='MINRES':                  if solver=='MINRES':      
1203                  if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1204                  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.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
# Line 1160  class SaddlePointProblem(object): Line 1272  class SaddlePointProblem(object):
1272         @param text: a text message         @param text: a text message
1273         @type text: C{str}         @type text: C{str}
1274         """         """
1275         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "#s: #s"%(str(self),text)
1276    
1277     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1278         """         """
Legend:
Removed from v.1488  
changed lines
  Added in v.1489
  ViewVC Help
Powered by ViewVC 1.1.26