/[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 1549 by artak, Wed Apr 30 03:40:15 2008 UTC revision 1550 by artak, Wed May 7 02:22:57 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,solver="GMRES",TOL=None):     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
# Line 493  class IterationHistory(object): Line 493  class IterationHistory(object):
493         if TOL==None:         if TOL==None:
494            TOL=self.tolerance            TOL=self.tolerance
495         self.history.append(norm_r)         self.history.append(norm_r)
496         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])
497         return self.history[-1]<=TOL * norm_b         return self.history[-1]<=TOL * norm_b
498    
499  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 581  def givapp(c,s,vin): Line 581  def givapp(c,s,vin):
581              vrot[i:i+2]=w1,w2              vrot[i:i+2]=w1,w2
582      return vrot      return vrot
583    
584    ##############################################
585  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
586    ################################################
587     m=iter_restart     m=iter_restart
588     iter=0     iter=0
589       xc=x
590     while True:     while True:
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        x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)        xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
593        if stopped: break        if stopped: break
594        iter+=iter_restart            iter+=iter_restart    
595     return x     return xc
596    
597  def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
598     iter=0     iter=0
# Line 599  def GMRESm(b, Aprod, Msolve, bilinearfor Line 602  def GMRESm(b, Aprod, Msolve, bilinearfor
602     norm_b=math.sqrt(r_dot_r)     norm_b=math.sqrt(r_dot_r)
603    
604     if x==None:     if x==None:
605        x=0        x=0*b
606     else:     else:
607        r=Msolve(b-Aprod(x))        r=Msolve(b-Aprod(x))
608        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
# Line 629  def GMRESm(b, Aprod, Msolve, bilinearfor Line 632  def GMRESm(b, Aprod, Msolve, bilinearfor
632  # Modified Gram-Schmidt  # Modified Gram-Schmidt
633      for j in range(iter+1):      for j in range(iter+1):
634        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j][iter]=bilinearform(v[j],v[iter+1])  
635        v[iter+1]+=(-1.)*h[j][iter]*v[j]        v[iter+1]-=h[j][iter]*v[j]
636                
637      h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))      h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
638      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1][iter]
# Line 640  def GMRESm(b, Aprod, Msolve, bilinearfor Line 643  def GMRESm(b, Aprod, Msolve, bilinearfor
643       for j in range(iter+1):       for j in range(iter+1):
644          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
645              h[j][iter]=h[j][iter]+hr #vhat              h[j][iter]=h[j][iter]+hr #vhat
646              v[iter+1] +=(-1.)*hr*v[j]              v[iter+1] -= hr*v[j]
647    
648       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
649       h[iter+1][iter]=v_norm2       h[iter+1][iter]=v_norm2
650    
651  #   watch out for happy breakdown  #   watch out for happy breakdown
652          if v_norm2 != 0:          if not v_norm2 == 0:
653           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1][iter]
654    
655  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
656      if iter > 0 :      if iter > 0 :
657          hhat=[]          hhat=numarray.zeros(iter+1,numarray.Float64)
658          for i in range(iter+1) : hhat.append(h[i][iter])          for i in range(iter+1) : hhat[i]=h[i][iter]
659          hhat=givapp(c[0:iter],s[0:iter],hhat);          hhat=givapp(c[0:iter],s[0:iter],hhat);
660              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i][iter]=hhat[i]
661    
# Line 692  def GMRESm(b, Aprod, Msolve, bilinearfor Line 695  def GMRESm(b, Aprod, Msolve, bilinearfor
695    
696     return x,stopped     return x,stopped
697    
698    
699  ######################################################  ######################################################
700  def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):  def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
701  ######################################################  ######################################################
# Line 734  def FDGMRES(f0, Aprod, Msolve, bilinearf Line 738  def FDGMRES(f0, Aprod, Msolve, bilinearf
738     r=b     r=b
739    
740     if x==None:     if x==None:
741        x=0        x=0*f0
742     else:     else:
743        r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0  
744                
# Line 827  def FDGMRES(f0, Aprod, Msolve, bilinearf Line 831  def FDGMRES(f0, Aprod, Msolve, bilinearf
831    
832     return x,stopped     return x,stopped
833    
834    #################################################
835  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
836    #################################################
837      #      #
838      #  minres solves the system of linear equations Ax = b      #  minres solves the system of linear equations Ax = b
839      #  where A is a symmetric matrix (possibly indefinite or singular)      #  where A is a symmetric matrix (possibly indefinite or singular)
# Line 858  def MINRES(b, Aprod, Msolve, bilinearfor Line 863  def MINRES(b, Aprod, Msolve, bilinearfor
863    
864      r1    = b      r1    = b
865      y = Msolve(b)      y = Msolve(b)
866      beta1 = bilinearform(b,y)      beta1 = bilinearform(y,b)
867    
868      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
869    
# Line 920  def MINRES(b, Aprod, Msolve, bilinearfor Line 925  def MINRES(b, Aprod, Msolve, bilinearfor
925            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
926    
927          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
928          y      = (- alfa/beta)*r2 + y          y      += (- alfa/beta)*r2
929          r1     = r2          r1     = r2
930          r2     = y          r2     = y
931          y = Msolve(r2)          y = Msolve(r2)
932          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
933          beta   = bilinearform(r2,y)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
934          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm,"negative norm."
935    
936          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
# Line 960  def MINRES(b, Aprod, Msolve, bilinearfor Line 965  def MINRES(b, Aprod, Msolve, bilinearfor
965          w1    = w2          w1    = w2
966          w2    = w          w2    = w
967          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
968          x     = x  +  phi*w          x     +=  phi*w
969    
970          # Go round again.          # Go round again.
971    
# Line 979  def MINRES(b, Aprod, Msolve, bilinearfor Line 984  def MINRES(b, Aprod, Msolve, bilinearfor
984          rnorm  = phibar          rnorm  = phibar
985    
986      return x      return x
       
 def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20):  
987    
988    ######################################    
989    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
990    #####################################
991      gamma=.9      gamma=.9
992      lmaxit=40      lmaxit=100
993      etamax=.5      etamax=.5
994    
995      n = 1 #len(x)      n = 1 #len(x)
# Line 1001  def NewtonGMRES(b, Aprod, Msolve, biline Line 1007  def NewtonGMRES(b, Aprod, Msolve, biline
1007      f0=-Msolve(r)      f0=-Msolve(r)
1008      fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)      fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1009      fnrmo=1      fnrmo=1
     atol=1.e-2  
     rtol=1.e-4  
1010      stop_tol=atol + rtol*fnrm      stop_tol=atol + rtol*fnrm
1011      #      #
1012      # main iteration loop      # main iteration loop
1013      #      #
1014      while not stoppingcriterium(fnrm,stop_tol,'NewtonGMRES',TOL=1.):      while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):
1015    
1016              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
1017          #          #
# Line 1021  def NewtonGMRES(b, Aprod, Msolve, biline Line 1025  def NewtonGMRES(b, Aprod, Msolve, biline
1025          #          #
1026              # compute the step using a GMRES(m) routine especially designed for this purpose              # compute the step using a GMRES(m) routine especially designed for this purpose
1027          #          #
1028              initer=0              initer=0
1029              while True:              while True:
1030                 xc,stopped=FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)                 xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1031                 if stopped: break                 if stopped: break
1032                 initer+=iter_restart                 initer+=iter_restart
         xold=x  
1033          x+=xc          x+=xc
1034          f0=-Msolve(Aprod(x))          f0=-Msolve(Aprod(x))
1035          fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)          fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
# Line 1265  class ArithmeticTuple(object): Line 1268  class ArithmeticTuple(object):
1268             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1269         return self         return self
1270    
1271       def __add__(self,other):
1272           """
1273           add other to self
1274    
1275           @param other: increment
1276           @type other: C{ArithmeticTuple}
1277           """
1278           if not isinstance(other,float):
1279            if len(self) != len(other):
1280               raise ValueError,"tuple length must match."
1281            for i in range(len(self)):
1282               self.__items[i]+=other[i]
1283           else:
1284            for i in range(len(self)):
1285               self.__items[i]+=other
1286    
1287           return self
1288    
1289       def __sub__(self,other):
1290           """
1291           subtract other from self
1292    
1293           @param other: increment
1294           @type other: C{ArithmeticTuple}
1295           """
1296           if len(self) != len(other):
1297               raise ValueError,"tuple length must match."
1298           for i in range(len(self)):
1299               self.__items[i]-=other[i]
1300           return self
1301    
1302       def __neg__(self):
1303           """
1304           negate
1305    
1306           """
1307           for i in range(len(self)):
1308               self.__items[i]=-self.__items[i]
1309           return self
1310    
1311    
1312  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1313        """        """
1314        This provides a framwork for solving homogeneous saddle point problem of the form        This provides a framwork for solving homogeneous saddle point problem of the form
# Line 1336  class HomogeneousSaddlePointProblem(obje Line 1380  class HomogeneousSaddlePointProblem(obje
1380    
1381        def __inner_p(self,p1,p2):        def __inner_p(self,p1,p2):
1382           return self.inner(p1,p2)           return self.inner(p1,p2)
1383          
1384          def __inner_a(self,a1,a2):
1385             return self.inner_a(a1,a2)
1386    
1387          def __inner_a1(self,a1,a2):
1388             return self.inner(a1[1],a2[1])
1389    
1390        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
1391            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
# Line 1389  class HomogeneousSaddlePointProblem(obje Line 1439  class HomogeneousSaddlePointProblem(obje
1439    
1440            if solver=='NewtonGMRES':                if solver=='NewtonGMRES':    
1441                  if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1442                  p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)                  p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,atol=0,rtol=self.getTolerance())
1443                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1444                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1445                  #       u=v+(u-v)                  #       u=v+(u-v)
# Line 1412  class HomogeneousSaddlePointProblem(obje Line 1462  class HomogeneousSaddlePointProblem(obje
1462                  #       u=v+(u-v)                  #       u=v+(u-v)
1463          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1464                                
1465              if solver=='GMRESC':      
1466                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1467                    p0=self.solve_prec(Bz)
1468               #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1469                   #p-=alfa
1470                    x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)
1471                    #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())
1472                    # solve Au=f-B^*p
1473                    #       A(u-v)=f-B^*p-Av
1474                    #       u=v+(u-v)
1475                p=x[1]
1476            #u=v+self.solve_A(v,p)      
1477            #p=x[1]
1478            #u=x[0]
1479    
1480                if solver=='PCG':                if solver=='PCG':
1481                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1482                  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)
# Line 1425  class HomogeneousSaddlePointProblem(obje Line 1490  class HomogeneousSaddlePointProblem(obje
1490            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1491    
1492        def __Msolve2(self,r):        def __Msolve2(self,r):
1493            return self.solve_prec(r)            return self.solve_prec(r*1)
1494    
1495          def __Mempty(self,r):
1496              return r
1497    
1498    
1499        def __Aprod(self,p):        def __Aprod(self,p):
# Line 1434  class HomogeneousSaddlePointProblem(obje Line 1502  class HomogeneousSaddlePointProblem(obje
1502            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1503            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1504    
1505          def __Anext(self,x):
1506              # return next v,p
1507              #solve Av =-B^*p as Av =f-Az-B^*p
1508    
1509          pc=x[1]
1510              v=self.solve_A(self.__z,-pc)
1511          p=self.solve_prec(self.B(v))
1512    
1513              return ArithmeticTuple(v,p)
1514    
1515    
1516        def __Aprod2(self,p):        def __Aprod2(self,p):
1517            # return BA^-1B*p            # return BA^-1B*p
1518            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
# Line 1441  class HomogeneousSaddlePointProblem(obje Line 1520  class HomogeneousSaddlePointProblem(obje
1520            return self.B(v)            return self.B(v)
1521    
1522        def __Aprod_Newton(self,p):        def __Aprod_Newton(self,p):
1523            # return BA^-1B*p            # return BA^-1B*p - Bz
1524            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
1525        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1526            return self.B(v-self.__z)            return self.B(v-self.__z)
1527    
1528          def __Aprod_Newton2(self,x):
1529              # return BA^-1B*p - Bz
1530              #solve Av =-B^*p as Av =f-Az-B^*p
1531              pc=x[1]
1532          v=self.solve_A(self.__z,-pc)
1533              p=self.solve_prec(self.B(v-self.__z))
1534              return ArithmeticTuple(v,p)
1535    
1536  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1537     """     """
1538     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 1483  class SaddlePointProblem(object): Line 1570  class SaddlePointProblem(object):
1570         @param text: a text message         @param text: a text message
1571         @type text: C{str}         @type text: C{str}
1572         """         """
1573         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "#s: #s"%(str(self),text)
1574    
1575     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1576         """         """
# Line 1633  class SaddlePointProblem(object): Line 1720  class SaddlePointProblem(object):
1720              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1721              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1722              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1723                self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
1724    
             self.trace("step %d: f = %e, f/u = %e, g = %e, g/p = %e, relaxation = %e."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))  
1725              if self.iter>1:              if self.iter>1:
1726                 dg2=g_new-g                 dg2=g_new-g
1727                 df2=f_new-f                 df2=f_new-f

Legend:
Removed from v.1549  
changed lines
  Added in v.1550

  ViewVC Help
Powered by ViewVC 1.1.26