/[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 1541 by ksteube, Wed Apr 30 02:05:29 2008 UTC revision 1552 by gross, Thu May 8 08:52:41 2008 UTC
# 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 1233  class ArithmeticTuple(object): Line 1236  class ArithmeticTuple(object):
1236         """         """
1237         out=[]         out=[]
1238         for i in range(len(self)):         for i in range(len(self)):
1239             out.append(self[i]/other)             out.append(self[i]*(1./other))
1240         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1241    
1242     def __rdiv__(self,other):     def __rdiv__(self,other):
# 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 1365  class HomogeneousSaddlePointProblem(obje Line 1415  class HomogeneousSaddlePointProblem(obje
1415                # which leads to BA^-1B^*p = BA^-1f                  # which leads to BA^-1B^*p = BA^-1f  
1416    
1417            # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)                  # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
   
             
1418            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
   
1419                Bz=self.B(self.__z)                Bz=self.B(self.__z)
1420                #                #
1421            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1422                #                #
               #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
1423                #                #
               #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)  
               #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)  
1424                #                #
1425                self.iter=0                self.iter=0
1426            if solver=='GMRES':                  if solver=='GMRES':      
# Line 1389  class HomogeneousSaddlePointProblem(obje Line 1433  class HomogeneousSaddlePointProblem(obje
1433    
1434            if solver=='NewtonGMRES':                if solver=='NewtonGMRES':    
1435                  if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1436                  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())
1437                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1438                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1439                  #       u=v+(u-v)                  #       u=v+(u-v)
# Line 1412  class HomogeneousSaddlePointProblem(obje Line 1456  class HomogeneousSaddlePointProblem(obje
1456                  #       u=v+(u-v)                  #       u=v+(u-v)
1457          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1458                                
1459              if solver=='GMRESC':      
1460                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1461                    p0=self.solve_prec(Bz)
1462               #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1463                   #p-=alfa
1464                    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)
1465                    #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())
1466                    # solve Au=f-B^*p
1467                    #       A(u-v)=f-B^*p-Av
1468                    #       u=v+(u-v)
1469                p=x[1]
1470            #u=v+self.solve_A(v,p)      
1471            #p=x[1]
1472            #u=x[0]
1473    
1474                if solver=='PCG':                if solver=='PCG':
1475                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1476                    #
1477                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1478                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1479                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1480                  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)
1481              u=r[0]                u=r[0]  
1482                    print "DIFF=",util.integrate(p)
1483    
1484                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1485    
# Line 1425  class HomogeneousSaddlePointProblem(obje Line 1489  class HomogeneousSaddlePointProblem(obje
1489            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1490    
1491        def __Msolve2(self,r):        def __Msolve2(self,r):
1492            return self.solve_prec(r)            return self.solve_prec(r*1)
1493    
1494          def __Mempty(self,r):
1495              return r
1496    
1497    
1498        def __Aprod(self,p):        def __Aprod(self,p):
1499            # return BA^-1B*p            # return BA^-1B*p
1500            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1501            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1502            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1503    
1504          def __Anext(self,x):
1505              # return next v,p
1506              #solve Av =-B^*p as Av =f-Az-B^*p
1507    
1508          pc=x[1]
1509              v=self.solve_A(self.__z,-pc)
1510          p=self.solve_prec(self.B(v))
1511    
1512              return ArithmeticTuple(v,p)
1513    
1514    
1515        def __Aprod2(self,p):        def __Aprod2(self,p):
1516            # return BA^-1B*p            # return BA^-1B*p
1517            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1518        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1519            return self.B(v)            return self.B(v)
1520    
1521        def __Aprod_Newton(self,p):        def __Aprod_Newton(self,p):
1522            # return BA^-1B*p            # return BA^-1B*p - Bz
1523            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
1524        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1525            return self.B(v-self.__z)            return self.B(v-self.__z)
1526    
1527          def __Aprod_Newton2(self,x):
1528              # return BA^-1B*p - Bz
1529              #solve Av =-B^*p as Av =f-Az-B^*p
1530              pc=x[1]
1531          v=self.solve_A(self.__z,-pc)
1532              p=self.solve_prec(self.B(v-self.__z))
1533              return ArithmeticTuple(v,p)
1534    
1535  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1536     """     """
1537     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 1483  class SaddlePointProblem(object): Line 1569  class SaddlePointProblem(object):
1569         @param text: a text message         @param text: a text message
1570         @type text: C{str}         @type text: C{str}
1571         """         """
1572         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "#s: #s"%(str(self),text)
1573    
1574     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1575         """         """
# Line 1633  class SaddlePointProblem(object): Line 1719  class SaddlePointProblem(object):
1719              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1720              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1721              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1722                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))
1723    
             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))  
1724              if self.iter>1:              if self.iter>1:
1725                 dg2=g_new-g                 dg2=g_new-g
1726                 df2=f_new-f                 df2=f_new-f

Legend:
Removed from v.1541  
changed lines
  Added in v.1552

  ViewVC Help
Powered by ViewVC 1.1.26