/[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 1527 by ksteube, Tue Apr 29 00:52:01 2008 UTC revision 1557 by artak, Mon May 19 04:44:27 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 615  def GMRESm(b, Aprod, Msolve, bilinearfor Line 618  def GMRESm(b, Aprod, Msolve, bilinearfor
618        
619     v.append(r/rho)     v.append(r/rho)
620     g[0]=rho     g[0]=rho
621    
622     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
623    
624      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
625    
       
626      p=Msolve(Aprod(v[iter]))      p=Msolve(Aprod(v[iter]))
627    
628      v.append(p)      v.append(p)
# Line 627  def GMRESm(b, Aprod, Msolve, bilinearfor Line 630  def GMRESm(b, Aprod, Msolve, bilinearfor
630      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
631    
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]
639    
   
640  # Reorthogonalize if needed  # Reorthogonalize if needed
641      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
642       for j in range(iter+1):       for j in range(iter+1):  
643          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
644              h[j][iter]=h[j][iter]+hr #vhat              h[j][iter]=h[j][iter]+hr
645              v[iter+1] +=(-1.)*hr*v[j]              v[iter+1] -= hr*v[j]
646    
647       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
648       h[iter+1][iter]=v_norm2       h[iter+1][iter]=v_norm2
649    
650  #   watch out for happy breakdown  #   watch out for happy breakdown
651          if v_norm2 != 0:          if not v_norm2 == 0:
652           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1][iter]
653    
654  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
655      if iter > 0 :      if iter > 0 :
656          hhat=[]          hhat=numarray.zeros(iter+1,numarray.Float64)
657          for i in range(iter+1) : hhat.append(h[i][iter])          for i in range(iter+1) : hhat[i]=h[i][iter]
658          hhat=givapp(c[0:iter],s[0:iter],hhat);          hhat=givapp(c[0:iter],s[0:iter],hhat);
659              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i][iter]=hhat[i]
660    
661      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
662    
663      if mu!=0 :      if mu!=0 :
664          c[iter]=h[iter][iter]/mu          c[iter]=h[iter][iter]/mu
665          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1][iter]/mu
# Line 665  def GMRESm(b, Aprod, Msolve, bilinearfor Line 668  def GMRESm(b, Aprod, Msolve, bilinearfor
668          g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])          g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
669    
670  # Update the residual norm  # Update the residual norm
671                  
672          rho=abs(g[iter+1])          rho=abs(g[iter+1])
673      iter+=1      iter+=1
674    
# Line 683  def GMRESm(b, Aprod, Msolve, bilinearfor Line 687  def GMRESm(b, Aprod, Msolve, bilinearfor
687       for i in range(iter-1):       for i in range(iter-1):
688      xhat += v[i]*y[i]      xhat += v[i]*y[i]
689     else : xhat=v[0]     else : xhat=v[0]
690        
691     x += xhat     x += xhat
692     if iter<iter_restart-1:     if iter<iter_restart-1:
693        stopped=True        stopped=True
# Line 692  def GMRESm(b, Aprod, Msolve, bilinearfor Line 696  def GMRESm(b, Aprod, Msolve, bilinearfor
696    
697     return x,stopped     return x,stopped
698    
699    
700  ######################################################  ######################################################
701  def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):  def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
702  ######################################################  ######################################################
# Line 734  def FDGMRES(f0, Aprod, Msolve, bilinearf Line 739  def FDGMRES(f0, Aprod, Msolve, bilinearf
739     r=b     r=b
740    
741     if x==None:     if x==None:
742        x=0        x=0*f0
743     else:     else:
744        r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0  
745                
# Line 827  def FDGMRES(f0, Aprod, Msolve, bilinearf Line 832  def FDGMRES(f0, Aprod, Msolve, bilinearf
832    
833     return x,stopped     return x,stopped
834    
835    #################################################
836  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
837    #################################################
838      #      #
839      #  minres solves the system of linear equations Ax = b      #  minres solves the system of linear equations Ax = b
840      #  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 864  def MINRES(b, Aprod, Msolve, bilinearfor
864    
865      r1    = b      r1    = b
866      y = Msolve(b)      y = Msolve(b)
867      beta1 = bilinearform(b,y)      beta1 = bilinearform(y,b)
868    
869      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
870    
# Line 920  def MINRES(b, Aprod, Msolve, bilinearfor Line 926  def MINRES(b, Aprod, Msolve, bilinearfor
926            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
927    
928          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
929          y      = (- alfa/beta)*r2 + y          y      += (- alfa/beta)*r2
930          r1     = r2          r1     = r2
931          r2     = y          r2     = y
932          y = Msolve(r2)          y = Msolve(r2)
933          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
934          beta   = bilinearform(r2,y)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
935          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm,"negative norm."
936    
937          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
# Line 960  def MINRES(b, Aprod, Msolve, bilinearfor Line 966  def MINRES(b, Aprod, Msolve, bilinearfor
966          w1    = w2          w1    = w2
967          w2    = w          w2    = w
968          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
969          x     = x  +  phi*w          x     +=  phi*w
970    
971          # Go round again.          # Go round again.
972    
# Line 979  def MINRES(b, Aprod, Msolve, bilinearfor Line 985  def MINRES(b, Aprod, Msolve, bilinearfor
985          rnorm  = phibar          rnorm  = phibar
986    
987      return x      return x
       
 def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20):  
988    
989    ######################################    
990    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
991    #####################################
992      gamma=.9      gamma=.9
993      lmaxit=40      lmaxit=100
994      etamax=.5      etamax=.5
995    
996      n = 1 #len(x)      n = 1 #len(x)
# Line 1001  def NewtonGMRES(b, Aprod, Msolve, biline Line 1008  def NewtonGMRES(b, Aprod, Msolve, biline
1008      f0=-Msolve(r)      f0=-Msolve(r)
1009      fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)      fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1010      fnrmo=1      fnrmo=1
     atol=1.e-2  
     rtol=1.e-4  
1011      stop_tol=atol + rtol*fnrm      stop_tol=atol + rtol*fnrm
1012      #      #
1013      # main iteration loop      # main iteration loop
1014      #      #
1015      while not stoppingcriterium(fnrm,stop_tol,'NewtonGMRES',TOL=1.):      while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):
1016    
1017              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
1018          #          #
# Line 1021  def NewtonGMRES(b, Aprod, Msolve, biline Line 1026  def NewtonGMRES(b, Aprod, Msolve, biline
1026          #          #
1027              # 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
1028          #          #
1029              initer=0              initer=0
1030              while True:              while True:
1031                 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)
1032                 if stopped: break                 if stopped: break
1033                 initer+=iter_restart                 initer+=iter_restart
         xold=x  
1034          x+=xc          x+=xc
1035          f0=-Msolve(Aprod(x))          f0=-Msolve(Aprod(x))
1036          fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)          fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
# Line 1201  class ArithmeticTuple(object): Line 1205  class ArithmeticTuple(object):
1205         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1206         """         """
1207         out=[]         out=[]
1208         for i in range(len(self)):         other=1.*other
1209           if isinstance(other,float):
1210        for i in range(len(self)):
1211             out.append(self[i]*other)             out.append(self[i]*other)
1212           else:
1213            for i in range(len(self)):
1214               out.append(self[i]*other[i])
1215         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1216    
1217     def __rmul__(self,other):     def __rmul__(self,other):
# Line 1215  class ArithmeticTuple(object): Line 1224  class ArithmeticTuple(object):
1224         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1225         """         """
1226         out=[]         out=[]
1227         for i in range(len(self)):         other=1.*other
1228           if isinstance(other,float):
1229        for i in range(len(self)):
1230             out.append(other*self[i])             out.append(other*self[i])
1231           else:
1232            for i in range(len(self)):
1233               out.append(other[i]*self[i])
1234         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1235    
1236  #########################  #########################
# Line 1232  class ArithmeticTuple(object): Line 1246  class ArithmeticTuple(object):
1246         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1247         """         """
1248         out=[]         out=[]
1249         for i in range(len(self)):         other=1.*other
1250             out.append(self[i]/other)         if isinstance(other,float):
1251        for i in range(len(self)):
1252               out.append(self[i]*(1./other))
1253           else:
1254            for i in range(len(self)):
1255               out.append(self[i]*(1./other[i]))
1256         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1257    
1258     def __rdiv__(self,other):     def __rdiv__(self,other):
# Line 1246  class ArithmeticTuple(object): Line 1265  class ArithmeticTuple(object):
1265         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1266         """         """
1267         out=[]         out=[]
1268         for i in range(len(self)):         other=1.*other
1269             out.append(other/self[i])         if isinstance(other,float):
1270            for i in range(len(self)):
1271               out.append(other*(1./self[i]))
1272           else:
1273            for i in range(len(self)):
1274               out.append(other[i]*(1./self[i]))
1275         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1276        
1277  ##########################################33  ##########################################33
# Line 1265  class ArithmeticTuple(object): Line 1289  class ArithmeticTuple(object):
1289             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1290         return self         return self
1291    
1292       def __add__(self,other):
1293           """
1294           add other to self
1295    
1296           @param other: increment
1297           @type other: C{ArithmeticTuple}
1298           """
1299    #      if not isinstance(other,float):
1300           if len(self) != len(other):
1301              raise ValueError,"tuple length must match."
1302           for i in range(len(self)):
1303              self.__items[i]+=other[i]
1304    #       else:
1305    #        for i in range(len(self)):
1306    #           self.__items[i]+=other
1307    
1308           return self
1309    
1310       def __sub__(self,other):
1311           """
1312           subtract other from self
1313    
1314           @param other: increment
1315           @type other: C{ArithmeticTuple}
1316           """
1317           if len(self) != len(other):
1318               raise ValueError,"tuple length must match."
1319           for i in range(len(self)):
1320               self.__items[i]-=other[i]
1321           return self
1322      
1323       def __isub__(self,other):
1324           """
1325           subtract other from self
1326    
1327           @param other: increment
1328           @type other: C{ArithmeticTuple}
1329           """
1330           if len(self) != len(other):
1331               raise ValueError,"tuple length must match."
1332           for i in range(len(self)):
1333               self.__items[i]-=other[i]
1334           return self
1335    
1336       def __neg__(self):
1337           """
1338           negate
1339    
1340           """
1341           for i in range(len(self)):
1342               self.__items[i]=-self.__items[i]
1343           return self
1344    
1345    
1346  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1347        """        """
1348        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 1414  class HomogeneousSaddlePointProblem(obje
1414    
1415        def __inner_p(self,p1,p2):        def __inner_p(self,p1,p2):
1416           return self.inner(p1,p2)           return self.inner(p1,p2)
1417          
1418          def __inner_a(self,a1,a2):
1419             return self.inner_a(a1,a2)
1420    
1421          def __inner_a1(self,a1,a2):
1422             return self.inner(a1[1],a2[1])
1423    
1424        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
1425            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
# Line 1365  class HomogeneousSaddlePointProblem(obje Line 1449  class HomogeneousSaddlePointProblem(obje
1449                # which leads to BA^-1B^*p = BA^-1f                  # which leads to BA^-1B^*p = BA^-1f  
1450    
1451            # 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)      
   
             
1452            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
1453                  Bz=self.B(self.__z)
               Bz=self.B(self.__z)  
1454                #                #
1455            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1456                #                #
               #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
1457                #                #
               #   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)  
1458                #                #
1459                self.iter=0                self.iter=0
1460            if solver=='GMRES':                  if solver=='GMRES':      
# Line 1389  class HomogeneousSaddlePointProblem(obje Line 1467  class HomogeneousSaddlePointProblem(obje
1467    
1468            if solver=='NewtonGMRES':                if solver=='NewtonGMRES':    
1469                  if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1470                  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())
1471                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1472                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1473                  #       u=v+(u-v)                  #       u=v+(u-v)
# Line 1412  class HomogeneousSaddlePointProblem(obje Line 1490  class HomogeneousSaddlePointProblem(obje
1490                  #       u=v+(u-v)                  #       u=v+(u-v)
1491          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1492                                
1493              if solver=='GMRESC':      
1494                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1495                    p0=self.solve_prec1(Bz)
1496                #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1497                    #p-=alfa
1498                    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)
1499                    #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())
1500    
1501                    # solve Au=f-B^*p
1502                    #       A(u-v)=f-B^*p-Av
1503                    #       u=v+(u-v)
1504                p=x[1]
1505            u=v+self.solve_A(v,p)      
1506            #p=x[1]
1507            #u=x[0]
1508    
1509                if solver=='PCG':                if solver=='PCG':
1510                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1511                    #
1512                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1513                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1514                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1515                  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)
1516              u=r[0]                u=r[0]  
1517                    print "DIFF=",util.integrate(p)
1518    
1519                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1520    
# Line 1425  class HomogeneousSaddlePointProblem(obje Line 1524  class HomogeneousSaddlePointProblem(obje
1524            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1525    
1526        def __Msolve2(self,r):        def __Msolve2(self,r):
1527            return self.solve_prec(r)            return self.solve_prec(r*1)
1528    
1529          def __Mempty(self,r):
1530              return r
1531    
1532    
1533        def __Aprod(self,p):        def __Aprod(self,p):
1534            # return BA^-1B*p            # return BA^-1B*p
1535            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1536            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1537            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1538    
1539          def __Anext(self,x):
1540              # return next v,p
1541              #solve Av =-B^*p as Av =f-Az-B^*p
1542    
1543          pc=x[1]
1544              v=self.solve_A(self.__z,-pc)
1545          p=self.solve_prec1(self.B(v))
1546    
1547              return ArithmeticTuple(v,p)
1548    
1549    
1550        def __Aprod2(self,p):        def __Aprod2(self,p):
1551            # return BA^-1B*p            # return BA^-1B*p
1552            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1553        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1554            return self.B(v)            return self.B(v)
1555    
1556        def __Aprod_Newton(self,p):        def __Aprod_Newton(self,p):
1557            # return BA^-1B*p            # return BA^-1B*p - Bz
1558            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
1559        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1560            return self.B(v-self.__z)            return self.B(v-self.__z)
1561    
1562          def __Aprod_Newton2(self,x):
1563              # return BA^-1B*p - Bz
1564              #solve Av =-B^*p as Av =f-Az-B^*p
1565              pc=x[1]
1566          v=self.solve_A(self.__z,-pc)
1567              p=self.solve_prec1(self.B(v-self.__z))
1568              return ArithmeticTuple(v,p)
1569    
1570  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1571     """     """
1572     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 1633  class SaddlePointProblem(object): Line 1754  class SaddlePointProblem(object):
1754              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1755              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1756              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1757                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))
1758    
             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))  
1759              if self.iter>1:              if self.iter>1:
1760                 dg2=g_new-g                 dg2=g_new-g
1761                 df2=f_new-f                 df2=f_new-f

Legend:
Removed from v.1527  
changed lines
  Added in v.1557

  ViewVC Help
Powered by ViewVC 1.1.26