/[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 2155 by gross, Wed Nov 26 08:13:00 2008 UTC revision 2156 by gross, Mon Dec 15 05:09:02 2008 UTC
# Line 442  class NegativeNorm(SolverSchemeException Line 442  class NegativeNorm(SolverSchemeException
442     """     """
443     pass     pass
444    
445  def PCG(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, initial_guess=True, verbose=False):  def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
446     """     """
447     Solver for     Solver for
448    
# Line 451  def PCG(b, Aprod, Msolve, bilinearform, Line 451  def PCG(b, Aprod, Msolve, bilinearform,
451     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
452     It uses the conjugate gradient method with preconditioner M providing an approximation of A.     It uses the conjugate gradient method with preconditioner M providing an approximation of A.
453    
454     The iteration is terminated if the C{stoppingcriterium} function return C{True}.     The iteration is terminated if  
455    
456       M{|r| <= atol+rtol*|r0|}
457    
458       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
459    
460       M{|r| =sqrt( bilinearform(Msolve(r),r))}
461    
462     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
463    
# Line 459  def PCG(b, Aprod, Msolve, bilinearform, Line 465  def PCG(b, Aprod, Msolve, bilinearform,
465     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
466     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst.
467    
468     @param b: the right hand side of the liner system. C{b} is altered.     @param r: initial residual M{r=b-Ax}. C{r} is altered.
469     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
470       @param x: an initial guess for the solution.
471       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
472     @param Aprod: returns the value Ax     @param Aprod: returns the value Ax
473     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
474         The returned object needs to be of the same type like argument C{b}.         The returned object needs to be of the same type like argument C{r}.
475     @param Msolve: solves Mx=r     @param Msolve: solves Mx=r
476     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}.     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{r}.
477        The returned object needs to be of the same type like argument C{x}.        The returned object needs to be of the same type like argument C{x}.
478     @param bilinearform: inner product C{<x,r>}     @param bilinearform: inner product C{<x,r>}
479     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
480         The returned value is a C{float}.         The returned value is a C{float}.
481     @param stoppingcriterium: function which returns True if a stopping criterium is meet.     @param atol: absolute tolerance
482         C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current     @type atol: non-negative C{float}
483         norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and     @param rtol: relative tolerance
484         the current solution approximation. C{stoppingcriterium} is called in each iteration step.     @type rtol: non-negative C{float}
    @type stoppingcriterium: function that returns C{True} or C{False}  
    @param x: an initial guess for the solution. If no C{x} is given 0*b is used.  
    @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)  
485     @param iter_max: maximum number of iteration steps.     @param iter_max: maximum number of iteration steps.
486     @type iter_max: C{int}     @type iter_max: C{int}
487     @return: the solution approximation and the corresponding residual     @return: the solution approximation and the corresponding residual
488     @rtype: C{tuple}     @rtype: C{tuple}
489     @warning: C{b} and C{x} are altered.     @warning: C{r} and C{x} are altered.
490     """     """
491     iter=0     iter=0
    if x==None:  
       x=0*b  
       r=b*1.  
       initial_guess=False  
    if initial_guess:  
       r=b+(-1)*Aprod(x)  
       x=x*1.  
    else:  
       x=x*0  
       r=b*1  
492     rhat=Msolve(r)     rhat=Msolve(r)
493     d = rhat     d = rhat
494     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
495     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm,"negative norm."
496     norm_b=math.sqrt(rhat_dot_r)     norm_r0=math.sqrt(rhat_dot_r)
497     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_b, atol+rtol*norm_b)     atol2=atol+rtol*norm_r0
498       if atol2<=0:
499          raise ValueError,"Non-positive tolarance."
500       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
501    
502       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
503      
504    
505     while not math.sqrt(rhat_dot_r) <= atol+rtol*norm_b:     while not math.sqrt(rhat_dot_r) <= atol2:
506         iter+=1         iter+=1
507         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
508    
# Line 795  def __FDGMRES(F0, defect, x0, atol, iter Line 796  def __FDGMRES(F0, defect, x0, atol, iter
796    
797     return xhat,stopped     return xhat,stopped
798    
799  ##############################################  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False):
800  def GMRES(b, Aprod, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, iter_restart=20, verbose=False):     """
801  ################################################     Solver for
802    
803       M{Ax=b}
804    
805       with a general operator A (more details required!).
806       It uses the generalized minimum residual method (GMRES).
807    
808       The iteration is terminated if  
809    
810       M{|r| <= atol+rtol*|r0|}
811    
812       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
813    
814       M{|r| =sqrt( bilinearform(r,r))}
815    
816       @param r: initial residual M{r=b-Ax}. C{r} is altered.
817       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
818       @param x: an initial guess for the solution.
819       @type x: same like C{r}
820       @param Aprod: returns the value Ax
821       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
822           The returned object needs to be of the same type like argument C{r}.
823       @param bilinearform: inner product C{<x,r>}
824       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
825           The returned value is a C{float}.
826       @param atol: absolute tolerance
827       @type atol: non-negative C{float}
828       @param rtol: relative tolerance
829       @type rtol: non-negative C{float}
830       @param iter_max: maximum number of iteration steps.
831       @type iter_max: C{int}
832       @param iter_restart: in order to ssave memory the orthogonalization process is terminated after C{iter_restart} steps and the
833                        iteration is restarted.
834       @type iter_restart: C{int}
835       @return: the solution approximation and the corresponding residual.
836       @rtype: C{tuple}
837       @warning: C{r} and C{x} are altered.
838       """
839     m=iter_restart     m=iter_restart
840       restarted=False
841     iter=0     iter=0
    if x == None:  
       xc=b*0  
    else:  
       xc=x*1  
842     if rtol>0:     if rtol>0:
843        r_dot_r = bilinearform(b, b)        r_dot_r = bilinearform(r, r)
844        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm,"negative norm."
845        atol2=atol+rtol*math.sqrt(r_dot_r)        atol2=atol+rtol*math.sqrt(r_dot_r)
846        if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)        if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
847     else:     else:
848        atol2=atol        atol2=atol
849        if verbose: print "GMRES: absolute tolerance = %e"%atol2        if verbose: print "GMRES: absolute tolerance = %e"%atol2
850       if atol2<=0:
851          raise ValueError,"Non-positive tolarance."
852        
853     while True:     while True:
854        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
855        xc,stopped=__GMRESm(b, Aprod, bilinearform, atol2, x=xc, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)        if restarted:
856             r2 = r-Aprod(x-x2)
857          else:
858             r2=1*r
859          x2=x*1.
860          x,stopped=__GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)
861        iter+=iter_restart            iter+=iter_restart    
862        if stopped: break        if stopped: break
863        if verbose: print "GMRES: restart."        if verbose: print "GMRES: restart."
864          restarted=True
865     if verbose: print "GMRES: tolerance has reached."     if verbose: print "GMRES: tolerance has reached."
866     return xc     return x
867    
868  def __GMRESm(b, Aprod, bilinearform, atol, x, iter_max=100, iter_restart=20, verbose=False):  def __GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):
869     iter=0     iter=0
870        
871     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
# Line 831  def __GMRESm(b, Aprod, bilinearform, ato Line 874  def __GMRESm(b, Aprod, bilinearform, ato
874     g=numarray.zeros(iter_restart+1,numarray.Float64)     g=numarray.zeros(iter_restart+1,numarray.Float64)
875     v=[]     v=[]
876    
    r=b-Aprod(x)  
877     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
878     if r_dot_r<0: raise NegativeNorm,"negative norm."     if r_dot_r<0: raise NegativeNorm,"negative norm."
879     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
# Line 913  def __GMRESm(b, Aprod, bilinearform, ato Line 955  def __GMRESm(b, Aprod, bilinearform, ato
955    
956     return x,stopped     return x,stopped
957    
958  #################################################  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
959  def MINRES(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):      """
960  #################################################      Solver for
961      #  
962      #  minres solves the system of linear equations Ax = b      M{Ax=b}
963      #  where A is a symmetric matrix (possibly indefinite or singular)  
964      #  and b is a given vector.      with a symmetric and positive definite operator A (more details required!).
965      #        It uses the minimum residual method (MINRES) with preconditioner M providing an approximation of A.
966      #  "A" may be a dense or sparse matrix (preferably sparse!)  
967      #  or the name of a function such that      The iteration is terminated if  
968      #               y = A(x)  
969      #  returns the product y = Ax for any given vector x.      M{|r| <= atol+rtol*|r0|}
970      #  
971      #  "M" defines a positive-definite preconditioner M = C C'.      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
972      #  "M" may be a dense or sparse matrix (preferably sparse!)  
973      #  or the name of a function such that      M{|r| =sqrt( bilinearform(Msolve(r),r))}
974      #  solves the system My = x for any given vector x.  
975      #      For details on the preconditioned conjugate gradient method see the book:
976      #  
977            Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
978        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
979        C. Romine, and H. van der Vorst.
980    
981        @param r: initial residual M{r=b-Ax}. C{r} is altered.
982        @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
983        @param x: an initial guess for the solution.
984        @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
985        @param Aprod: returns the value Ax
986        @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
987            The returned object needs to be of the same type like argument C{r}.
988        @param Msolve: solves Mx=r
989        @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{r}.
990            The returned object needs to be of the same type like argument C{x}.
991        @param bilinearform: inner product C{<x,r>}
992        @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
993           The returned value is a C{float}.
994        @param atol: absolute tolerance
995        @type atol: non-negative C{float}
996        @param rtol: relative tolerance
997        @type rtol: non-negative C{float}
998        @param iter_max: maximum number of iteration steps.
999        @type iter_max: C{int}
1000        @return: the solution approximation and the corresponding residual
1001        @rtype: C{tuple}
1002        @warning: C{r} and C{x} are altered.
1003        """
1004      #------------------------------------------------------------------      #------------------------------------------------------------------
1005      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
1006      # y  =  beta1 P' v1,  where  P = C**(-1).      # y  =  beta1 P' v1,  where  P = C**(-1).
1007      # v is really P' v1.      # v is really P' v1.
1008      #------------------------------------------------------------------      #------------------------------------------------------------------
1009      if x==None:      r1    = r
1010        x=0*b      y = Msolve(r)
1011      else:      beta1 = bilinearform(y,r)
       b += (-1)*Aprod(x)  
   
     r1    = b  
     y = Msolve(b)  
     beta1 = bilinearform(y,b)  
1012    
1013      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
1014    
1015      #  If b = 0 exactly, stop with x = 0.      #  If r = 0 exactly, stop with x
1016      if beta1==0: return x*0.      if beta1==0: return x
1017    
1018      if beta1> 0:      if beta1> 0: beta1  = math.sqrt(beta1)      
       beta1  = math.sqrt(beta1)        
1019    
1020      #------------------------------------------------------------------      #------------------------------------------------------------------
1021      # Initialize quantities.      # Initialize quantities.
# Line 973  def MINRES(b, Aprod, Msolve, bilinearfor Line 1035  def MINRES(b, Aprod, Msolve, bilinearfor
1035      ynorm2 = 0      ynorm2 = 0
1036      cs     = -1      cs     = -1
1037      sn     = 0      sn     = 0
1038      w      = b*0.      w      = r*0.
1039      w2     = b*0.      w2     = r*0.
1040      r2     = r1      r2     = r1
1041      eps    = 0.0001      eps    = 0.0001
1042    
# Line 1067  def MINRES(b, Aprod, Msolve, bilinearfor Line 1129  def MINRES(b, Aprod, Msolve, bilinearfor
1129    
1130      return x      return x
1131    
1132  def TFQMR(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):  def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1133      """
1134  # TFQMR solver for linear systems    Solver for
1135  #  
1136  #    M{Ax=b}
1137  # initialization  
1138  #    with a general operator A (more details required!).
1139    errtol = math.sqrt(bilinearform(b,b))    It uses the generalized minimum residual method (GMRES).
1140    norm_b=errtol  
1141    kmax  = iter_max    The iteration is terminated if  
1142    error = []  
1143      M{|r| <= atol+rtol*|r0|}
1144    if math.sqrt(bilinearform(x,x)) != 0.0:  
1145      r = b - Aprod(x)    where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1146    else:  
1147      r = b    M{|r| =sqrt( bilinearform(r,r))}
1148    
1149    r=Msolve(r)    @param r: initial residual M{r=b-Ax}. C{r} is altered.
1150      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1151      @param x: an initial guess for the solution.
1152      @type x: same like C{r}
1153      @param Aprod: returns the value Ax
1154      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
1155          The returned object needs to be of the same type like argument C{r}.
1156      @param bilinearform: inner product C{<x,r>}
1157      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
1158          The returned value is a C{float}.
1159      @param atol: absolute tolerance
1160      @type atol: non-negative C{float}
1161      @param rtol: relative tolerance
1162      @type rtol: non-negative C{float}
1163      @param iter_max: maximum number of iteration steps.
1164      @type iter_max: C{int}
1165      @rtype: C{tuple}
1166      @warning: C{r} and C{x} are altered.
1167      """
1168    u1=0    u1=0
1169    u2=0    u2=0
1170    y1=0    y1=0
# Line 1095  def TFQMR(b, Aprod, Msolve, bilinearform Line 1174  def TFQMR(b, Aprod, Msolve, bilinearform
1174    y1 = r    y1 = r
1175    iter = 0    iter = 0
1176    d = 0    d = 0
1177        v = Aprod(y1)
   v = Msolve(Aprod(y1))  
1178    u1 = v    u1 = v
1179        
1180    theta = 0.0;    theta = 0.0;
1181    eta = 0.0;    eta = 0.0;
1182    tau = math.sqrt(bilinearform(r,r))    rho=bilinearform(r,r)
1183    error = [ error, tau ]    if rho < 0: raise NegativeNorm,"negative norm."
1184    rho = tau * tau    tau = math.sqrt(rho)
1185  #    norm_r0=tau
1186  #  TFQMR iteration    while tau>atol+rtol*norm_r0:
 #  
 #  while ( iter < kmax-1 ):  
     
   while not tau<=atol+rtol*norm_b:  
1187      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
1188    
1189      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1190        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
     if ( sigma == 0.0 ):  
       raise 'TFQMR breakdown, sigma=0'  
       
   
1191      alpha = rho / sigma      alpha = rho / sigma
1192    
1193      for j in range(2):      for j in range(2):
# Line 1126  def TFQMR(b, Aprod, Msolve, bilinearform Line 1196  def TFQMR(b, Aprod, Msolve, bilinearform
1196  #  #
1197        if ( j == 1 ):        if ( j == 1 ):
1198          y2 = y1 - alpha * v          y2 = y1 - alpha * v
1199          u2 = Msolve(Aprod(y2))          u2 = Aprod(y2)
1200    
1201        m = 2 * (iter+1) - 2 + (j+1)        m = 2 * (iter+1) - 2 + (j+1)
1202        if j==0:        if j==0:
# Line 1144  def TFQMR(b, Aprod, Msolve, bilinearform Line 1214  def TFQMR(b, Aprod, Msolve, bilinearform
1214  #  #
1215  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1216  #  #
1217       # if ( tau * math.sqrt ( m + 1 ) <= errtol ):      if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
      #   error = [ error, tau ]  
      #   total_iters = iter  
      #   break  
         
   
     if ( rho == 0.0 ):  
       raise 'TFQMR breakdown, rho=0'  
       
1218    
1219      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1220      beta = rhon / rho;      beta = rhon / rho;
1221      rho = rhon;      rho = rhon;
1222      y1 = w + beta * y2;      y1 = w + beta * y2;
1223      u1 = Msolve(Aprod(y1))      u1 = Aprod(y1)
1224      v = u1 + beta * ( u2 + beta * v )      v = u1 + beta * ( u2 + beta * v )
     error = [ error, tau ]  
     total_iters = iter  
1225            
1226      iter = iter + 1      iter += 1
1227    
1228    return x    return x
1229    
# Line 1363  class HomogeneousSaddlePointProblem(obje Line 1423  class HomogeneousSaddlePointProblem(obje
1423        def __init__(self,**kwargs):        def __init__(self,**kwargs):
1424          self.setTolerance()          self.setTolerance()
1425          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1426          self.setSubToleranceReductionFactor()          self.setSubProblemTolerance()
1427    
1428        #=============================================================        #=============================================================
1429        def initialize(self):        def initialize(self):
# Line 1494  class HomogeneousSaddlePointProblem(obje Line 1554  class HomogeneousSaddlePointProblem(obje
1554           self.show_details=show_details and self.verbose           self.show_details=show_details and self.verbose
1555           #=================================================================================           #=================================================================================
1556           # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)           # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)
          self.setSubProblemTolerance(self.getTolerance(), apply_reduction=True)  
1557           self.__z=v+self.solve_A(v,p*0)           self.__z=v+self.solve_A(v,p*0)
          Bz=self.B(self.__z)  
          p0=self.solve_prec(Bz)  
1558           # tolerances:           # tolerances:
1559           rtol=self.getTolerance()           rtol=self.getTolerance()
1560           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1561           f0=self.norm(self.__z,p0)           if useUzawa:
1562           if not f0>0:              p0=self.solve_prec(self.B(self.__z))
1563              return self.__z, p*0              f0=self.norm(self.__z,p0)
1564             else:
1565                f0=util.sqrt(self.inner_v(self.__z,self.__z))
1566             if not f0>0: return self.__z, p*0
1567           ATOL=rtol*f0+atol           ATOL=rtol*f0+atol
1568           if not ATOL>0:           if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
             raise ValueError,"overall absolute tolerance needs to be positive."  
1569           if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)           if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)
1570           # initialization           # initialization
1571           self.iter=0           self.iter=0
# Line 1514  class HomogeneousSaddlePointProblem(obje Line 1573  class HomogeneousSaddlePointProblem(obje
1573           converged=False           converged=False
1574           # initial guess:           # initial guess:
1575           q=p*1           q=p*1
1576           if not useUzawa:           u=v*1
            u=v  
1577           while not converged :           while not converged :
1578              if useUzawa:              if useUzawa:
                self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)  
1579                 # assume p is known: then v=z-A^-1B^*p                 # assume p is known: then v=z-A^-1B^*p
1580                 #                       #      
1581                 # which leads to BA^-1B^*p = Bz                 # which leads to BA^-1B^*p = Bz
1582                 #                 #
1583                 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bu                 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1584                   # we use the (v,Bv) to represent the residual
1585                 #                 #
1586                 # note that q=BA^-1B^*p = B (A^-1B^*p)=Bw  with Aw=B^*p = f - Az - B^*(-p)                 # the norm of the right hand side Bv = f0
1587                 #                 #
1588                   #                  and the initial residual
1589                   #          
1590                   #     r=Bz-BA^-1B^*q = B(z-A^{-1}B^*q)=Bw with A(w-z)=Az-Az-B^*q = f -A*0 - B^{*}q
1591                   #
1592                   w=self.solve_A(self.__z,q)+self.__z
1593                 if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())                 if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())
1594                 q,r=PCG(ArithmeticTuple(self.__z,Bz),self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=q, verbose=self.verbose)                 q,r=PCG(ArithmeticTuple(w,self.B(w)),self.__Aprod_PCG,q,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1595             u=r[0]               u=r[0]  
1596             Bu=r[1]               Bu=r[1]
1597              else:              else:
                self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)  
1598                 #                 #
1599                 #  with v=dv+z                 #  with v=dv+z
1600                 #                 #
1601                 #   A dv + B p = A (v-z) + Bp = 0                 #   A v + B p = f
1602                 #   B dv       =  - Bz                 #   B v       = 0
1603                 #                 #
1604                 # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1})  -S^{-1}]] to the problem                 # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1})  -S^{-1}]]
                # the right hand side is then [0, S^{-1} Bz ]  
1605                 #                 #
1606                 #                           w=self.solve_A(u,q)
1607                 if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())                 if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())
1608                 x=GMRES(ArithmeticTuple(0*self.__z,p0),self.__Aprod_GMRES,self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, x=ArithmeticTuple(u-self.__z,q),iter_restart=iter_restart, verbose=self.verbose)                 x=GMRES(ArithmeticTuple(w,self.solve_prec(self.B(u+w))),self.__Aprod_GMRES, ArithmeticTuple(u,q), \
1609                 self.setSubProblemTolerance(ATOL/f0, apply_reduction=False)                           self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1610             u=self.__z+x[0]             u=x[0]
1611                 q=x[1]                 q=x[1]
1612             Bu=self.B(u)             Bu=self.B(u)
1613              # now we check if |Bu| ~ 0 or more precise |Bu|  <= rtol * |v|              # now we check if |Bu| ~ 0 or more precise |Bu|_p  <= rtol * |v|_v
1614              nu=self.inner_v(u,u)              nu=self.inner_v(u,u)
             nz=self.inner_v(self.__z,self.__z)  
1615              p2=self.solve_prec(Bu)              p2=self.solve_prec(Bu)
1616              nBu=self.inner_p(p2,p2)              nBu=self.inner_p(p2,p2)
1617              if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."              if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."
# Line 1582  class HomogeneousSaddlePointProblem(obje Line 1642  class HomogeneousSaddlePointProblem(obje
1642           if tolerance<0:           if tolerance<0:
1643               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
1644           self.__rtol=tolerance           self.__rtol=tolerance
1645             self.setSubProblemTolerance()
1646    
1647        def getTolerance(self):        def getTolerance(self):
1648           """           """
1649           returns the relative tolerance           returns the relative tolerance
# Line 1608  class HomogeneousSaddlePointProblem(obje Line 1670  class HomogeneousSaddlePointProblem(obje
1670           @rtype: C{float}           @rtype: C{float}
1671           """           """
1672           return self.__atol           return self.__atol
       def setSubToleranceReductionFactor(self,reduction=None):  
          """  
          sets the reduction factor for the tolerance used to solve the subproblem in each iteration step.  
          if set to C{None} the square of the problem tolerance is used.  
       
          @param reduction: reduction factor.  
          @type reduction: positive C{float} less or equal one or C{None}  
          """  
          if not reduction == None:  
               if reduction<=0 or reduction>1:  
                   raise ValueError,"reduction factor (=%e) must be positive and less or equal one."%reduction  
          self.__reduction=reduction  
       def getSubToleranceReductionFactor(self):  
          """  
          returns the subproblem reduction factor  
   
          @return: subproblem reduction factor  
          @rtype: C{float} or C{None}  
          """  
          return self.__reduction  
1673    
1674        def setSubProblemTolerance(self,rtol,apply_reduction=False):        def setSubProblemTolerance(self,rtol=None):
1675           """           """
1676           sets the relative tolerance to solve the subproblem(s).           sets the relative tolerance to solve the subproblem(s).
1677    
1678           @param rtol: relative tolerence           @param rtol: relative tolerence
1679           @type rtol: positive C{float}           @type rtol: positive C{float}
          @param apply_reduction: if set the the reduction factor is applied.  
          @type apply_reduction: C{bool}  
1680           """           """
1681             if rtol == None:
1682                  rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1683           if rtol<=0:           if rtol<=0:
1684               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
1685           if apply_reduction:           self.__sub_tol=rtol
              if self.__reduction == None:  
                  self.__sub_tol=max(rtol**2,util.EPSILON*10)  
              else:  
                  self.__sub_tol=max(self.__reduction*rtol,util.EPSILON*10)  
          else:  
              self.__sub_tol=rtol  
1686        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1687           """           """
1688           returns the subproblem reduction factor           returns the subproblem reduction factor

Legend:
Removed from v.2155  
changed lines
  Added in v.2156

  ViewVC Help
Powered by ViewVC 1.1.26