/[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 1956 by gross, Mon Nov 3 05:08:42 2008 UTC revision 2156 by gross, Mon Dec 15 05:09:02 2008 UTC
# Line 426  class MaxIterReached(SolverSchemeExcepti Line 426  class MaxIterReached(SolverSchemeExcepti
426     maxium number of iteration steps is reached.     maxium number of iteration steps is reached.
427     """     """
428     pass     pass
429    class CorrectionFailed(SolverSchemeException):
430       """
431       no convergence has been achieved in the solution correction scheme.
432       """
433       pass
434  class IterationBreakDown(SolverSchemeException):  class IterationBreakDown(SolverSchemeException):
435     """     """
436     iteration scheme econouters an incurable breakdown.     iteration scheme econouters an incurable breakdown.
# Line 437  class NegativeNorm(SolverSchemeException Line 442  class NegativeNorm(SolverSchemeException
442     """     """
443     pass     pass
444    
445  class IterationHistory(object):  def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
    """  
    The IterationHistory class is used to define a stopping criterium. It keeps track of the  
    residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by  
    a given tolerance.  
    """  
    def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):  
       """  
       Initialization  
   
       @param tolerance: tolerance  
       @type tolerance: positive C{float}  
       @param verbose: switches on the printing out some information  
       @type verbose: C{bool}  
       """  
       if not tolerance>0.:  
           raise ValueError,"tolerance needs to be positive."  
       self.tolerance=tolerance  
       self.verbose=verbose  
       self.history=[]  
    def stoppingcriterium(self,norm_r,r,x):  
        """  
        returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.  
   
         
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param r: current residual (not used)  
        @param x: current solution approximation (not used)  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
   
        """  
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=self.tolerance * self.history[0]  
   
    def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):  
        """  
        returns True if the C{norm_r} is C{tolerance}*C{norm_b}  
   
         
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param norm_b: norm of right hand side  
        @type norm_b: non-negative C{float}  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
   
        """  
        if TOL==None:  
           TOL=self.tolerance  
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=TOL * norm_b  
   
 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  
446     """     """
447     Solver for     Solver for
448    
# Line 502  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 510  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 Aprod: returns the value Ax     @param x: an initial guess for the solution.
    @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.  
    @param Msolve: solves Mx=r  
    @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same  
 type like argument C{x}.  
    @param bilinearform: inner product C{<x,r>}  
    @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.  
    @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.  
    @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.  
471     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
472       @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}.
474           The returned object needs to be of the same type like argument C{r}.
475       @param Msolve: solves Mx=r
476       @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}.
478       @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.
480           The returned value is a C{float}.
481       @param atol: absolute tolerance
482       @type atol: non-negative C{float}
483       @param rtol: relative tolerance
484       @type rtol: non-negative C{float}
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  
    else:  
       b += (-1)*Aprod(x)  
    r=b  
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_r0=math.sqrt(rhat_dot_r)
497       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     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
503      
504    
505       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 557  type like argument C{x}. Line 519  type like argument C{x}.
519    
520         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
521         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
522           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
523       if verbose: print "PCG: tolerance reached after %s steps."%iter
524     return x,r     return x,r
525    
526  class Defect(object):  class Defect(object):
# Line 771  def __FDGMRES(F0, defect, x0, atol, iter Line 734  def __FDGMRES(F0, defect, x0, atol, iter
734    
735          # Modified Gram-Schmidt          # Modified Gram-Schmidt
736      for j in range(iter+1):      for j in range(iter+1):
737           h[j][iter]=defect.bilinearform(v[j],v[iter+1])             h[j,iter]=defect.bilinearform(v[j],v[iter+1])  
738           v[iter+1]-=h[j][iter]*v[j]           v[iter+1]-=h[j,iter]*v[j]
739                
740      h[iter+1][iter]=defect.norm(v[iter+1])      h[iter+1,iter]=defect.norm(v[iter+1])
741      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1,iter]
742    
743          # Reorthogonalize if needed          # Reorthogonalize if needed
744      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)
745          for j in range(iter+1):            for j in range(iter+1):  
746             hr=defect.bilinearform(v[j],v[iter+1])             hr=defect.bilinearform(v[j],v[iter+1])
747                 h[j][iter]=h[j][iter]+hr                 h[j,iter]=h[j,iter]+hr
748                 v[iter+1] -= hr*v[j]                 v[iter+1] -= hr*v[j]
749    
750          v_norm2=defect.norm(v[iter+1])          v_norm2=defect.norm(v[iter+1])
751          h[iter+1][iter]=v_norm2          h[iter+1,iter]=v_norm2
752          #   watch out for happy breakdown          #   watch out for happy breakdown
753          if not v_norm2 == 0:          if not v_norm2 == 0:
754                  v[iter+1]=v[iter+1]/h[iter+1][iter]                  v[iter+1]=v[iter+1]/h[iter+1,iter]
755    
756          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
757      if iter > 0 :      if iter > 0 :
758          hhat=numarray.zeros(iter+1,numarray.Float64)          hhat=numarray.zeros(iter+1,numarray.Float64)
759          for i in range(iter+1) : hhat[i]=h[i][iter]          for i in range(iter+1) : hhat[i]=h[i,iter]
760          hhat=__givapp(c[0:iter],s[0:iter],hhat);          hhat=__givapp(c[0:iter],s[0:iter],hhat);
761              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
762    
763      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])
764    
765      if mu!=0 :      if mu!=0 :
766          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
767          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
768          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
769          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
770          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])
771    
772          # Update the residual norm          # Update the residual norm
# Line 814  def __FDGMRES(F0, defect, x0, atol, iter Line 777  def __FDGMRES(F0, defect, x0, atol, iter
777     # It's time to compute x and leave.             # It's time to compute x and leave.        
778     if iter > 0 :     if iter > 0 :
779       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)    
780       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
781       if iter > 1 :         if iter > 1 :  
782          i=iter-2            i=iter-2  
783          while i>=0 :          while i>=0 :
784            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
785            i=i-1            i=i-1
786       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
787       for i in range(iter-1):       for i in range(iter-1):
# Line 833  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, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):     """
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
842     xc=x     if rtol>0:
843          r_dot_r = bilinearform(r, r)
844          if r_dot_r<0: raise NegativeNorm,"negative norm."
845          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)
847       else:
848          atol2=atol
849          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*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)        if restarted:
856        if stopped: break           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     return xc        if stopped: break
863          if verbose: print "GMRES: restart."
864          restarted=True
865       if verbose: print "GMRES: tolerance has reached."
866       return x
867    
868  def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def __GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):
869     iter=0     iter=0
    r=Msolve(b)  
    r_dot_r = bilinearform(r, r)  
    if r_dot_r<0: raise NegativeNorm,"negative norm."  
    norm_b=math.sqrt(r_dot_r)  
   
    if x==None:  
       x=0*b  
    else:  
       r=Msolve(b-Aprod(x))  
       r_dot_r = bilinearform(r, r)  
       if r_dot_r<0: raise NegativeNorm,"negative norm."  
870        
871     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
872     c=numarray.zeros(iter_restart,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
873     s=numarray.zeros(iter_restart,numarray.Float64)     s=numarray.zeros(iter_restart,numarray.Float64)
874     g=numarray.zeros(iter_restart,numarray.Float64)     g=numarray.zeros(iter_restart+1,numarray.Float64)
875     v=[]     v=[]
876    
877       r_dot_r = bilinearform(r, r)
878       if r_dot_r<0: raise NegativeNorm,"negative norm."
879     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
880        
881     v.append(r/rho)     v.append(r/rho)
882     g[0]=rho     g[0]=rho
883    
884     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
885       while not (rho<=atol or iter==iter_restart):
886    
887      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
888    
889      p=Msolve(Aprod(v[iter]))      p=Aprod(v[iter])
   
890      v.append(p)      v.append(p)
891    
892      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
893    
894  # Modified Gram-Schmidt  # Modified Gram-Schmidt
895      for j in range(iter+1):      for j in range(iter+1):
896        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j,iter]=bilinearform(v[j],v[iter+1])  
897        v[iter+1]-=h[j][iter]*v[j]        v[iter+1]-=h[j,iter]*v[j]
898                
899      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]))
900      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1,iter]
901    
902  # Reorthogonalize if needed  # Reorthogonalize if needed
903      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)
904       for j in range(iter+1):         for j in range(iter+1):  
905          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
906              h[j][iter]=h[j][iter]+hr              h[j,iter]=h[j,iter]+hr
907              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
908    
909       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
910       h[iter+1][iter]=v_norm2       h[iter+1,iter]=v_norm2
911    
912  #   watch out for happy breakdown  #   watch out for happy breakdown
913          if not v_norm2 == 0:          if not v_norm2 == 0:
914           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
915    
916  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
917      if iter > 0 :      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
918          hhat=numarray.zeros(iter+1,numarray.Float64)      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
         for i in range(iter+1) : hhat[i]=h[i][iter]  
         hhat=__givapp(c[0:iter],s[0:iter],hhat);  
             for i in range(iter+1) : h[i][iter]=hhat[i]  
   
     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])  
919    
920      if mu!=0 :      if mu!=0 :
921          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
922          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
923          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
924          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
925          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])
   
926  # Update the residual norm  # Update the residual norm
927                                
928          rho=abs(g[iter+1])          rho=abs(g[iter+1])
929            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
930      iter+=1      iter+=1
931    
932  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
933  # It's time to compute x and leave.          # It's time to compute x and leave.        
934    
935       if verbose: print "GMRES: iteration stopped after %s step."%iter
936     if iter > 0 :     if iter > 0 :
937       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)    
938       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
939       if iter > 1 :         if iter > 1 :  
940          i=iter-2            i=iter-2  
941          while i>=0 :          while i>=0 :
942            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
943            i=i-1            i=i-1
944       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
945       for i in range(iter-1):       for i in range(iter-1):
946      xhat += v[i]*y[i]      xhat += v[i]*y[i]
947     else : xhat=v[0]     else:
948         xhat=v[0] * 0
949    
950     x += xhat     x += xhat
951     if iter<iter_restart-1:     if iter<iter_restart-1:
# Line 948  def __GMRESm(b, Aprod, Msolve, bilinearf Line 955  def __GMRESm(b, Aprod, Msolve, bilinearf
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, stoppingcriterium, 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 1008  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    
1043      #---------------------------------------------------------------------      #---------------------------------------------------------------------
1044      # Main iteration loop.      # Main iteration loop.
1045      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1046      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1047    
1048      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
1049          iter    = iter  +  1          iter    = iter  +  1
# Line 1102  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, stoppingcriterium, 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 1130  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    m=1    norm_r0=tau
1186  #    while tau>atol+rtol*norm_r0:
 #  TFQMR iteration  
 #  
 #  while ( iter < kmax-1 ):  
     
   while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):  
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 1162  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 1180  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 1394  class HomogeneousSaddlePointProblem(obje Line 1418  class HomogeneousSaddlePointProblem(obje
1418               Bv    =0               Bv    =0
1419    
1420        for the unknowns v and p and given operators A and B and given right hand side f.        for the unknowns v and p and given operators A and B and given right hand side f.
1421        B^* is the adjoint operator of B is the given inner product.        B^* is the adjoint operator of B.
   
1422        """        """
1423        def __init__(self,**kwargs):        def __init__(self,**kwargs):
1424          self.setTolerance()          self.setTolerance()
1425          self.setToleranceReductionFactor()          self.setAbsoluteTolerance()
1426            self.setSubProblemTolerance()
1427    
1428          #=============================================================
1429        def initialize(self):        def initialize(self):
1430          """          """
1431          initialize the problem (overwrite)          initialize the problem (overwrite)
1432          """          """
1433          pass          pass
1434    
1435        def B(self,v):        def B(self,v):
1436            """
1437            returns Bv (overwrite)
1438            @rtype: equal to the type of p
1439    
1440            @note: boundary conditions on p should be zero!
1441            """
1442            raise NotImplementedError,"no operator B implemented."
1443    
1444          def inner_pBv(self,p,Bv):
1445           """           """
1446           returns Bv (overwrite)           returns inner product of element p and Bv  (overwrite)
1447           @rtype: equal to the type of p          
1448             @type p: equal to the type of p
1449             @type Bv: equal to the type of result of operator B
1450             @rtype: C{float}
1451    
1452           @note: boundary conditions on p should be zero!           @rtype: equal to the type of p
1453           """           """
1454           pass           raise NotImplementedError,"no inner product for p implemented."
1455    
1456        def inner(self,p0,p1):        def inner_p(self,p0,p1):
1457           """           """
1458           returns inner product of two element p0 and p1  (overwrite)           returns inner product of element p0 and p1  (overwrite)
1459                    
1460           @type p0: equal to the type of p           @type p0: equal to the type of p
1461           @type p1: equal to the type of p           @type p1: equal to the type of p
1462             @rtype: equal to the type of p
1463             """
1464             raise NotImplementedError,"no inner product for p implemented."
1465    
1466          def inner_v(self,v0,v1):
1467             """
1468             returns inner product of two element v0 and v1  (overwrite)
1469            
1470             @type v0: equal to the type of v
1471             @type v1: equal to the type of v
1472           @rtype: C{float}           @rtype: C{float}
1473    
1474           @rtype: equal to the type of p           @rtype: equal to the type of v
1475           """           """
1476             raise NotImplementedError,"no inner product for v implemented."
1477           pass           pass
1478    
1479        def solve_A(self,u,p):        def solve_A(self,u,p):
1480           """           """
1481           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           solves Av=f-Au-B^*p with accuracy self.getSubProblemTolerance() (overwrite)
1482    
1483           @rtype: equal to the type of v           @rtype: equal to the type of v
1484           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
1485           """           """
1486           pass           raise NotImplementedError,"no operator A implemented."
1487    
1488        def solve_prec(self,p):        def solve_prec(self,p):
1489           """           """
1490           provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)           provides a preconditioner for BA^{-1}B^* with accuracy self.getSubProblemTolerance() (overwrite)
1491    
1492           @rtype: equal to the type of p           @rtype: equal to the type of p
1493             @note: boundary conditions on p should be zero!
1494           """           """
1495           pass           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1496          #=============================================================
1497        def stoppingcriterium(self,Bv,v,p):        def __Aprod_PCG(self,p):
1498           """            # return (v,Bv) with v=A^-1B*p
1499           returns a True if iteration is terminated. (overwrite)            #solve Av =B^*p as Av =f-Az-B^*(-p)
1500              v=self.solve_A(self.__z,-p)
1501           @rtype: C{bool}            return ArithmeticTuple(v, self.B(v))
          """  
          pass  
               
       def __inner(self,p,r):  
          return self.inner(p,r[1])  
   
       def __inner_p(self,p1,p2):  
          return self.inner(p1,p2)  
         
       def __inner_a(self,a1,a2):  
          return self.inner_a(a1,a2)  
   
       def __inner_a1(self,a1,a2):  
          return self.inner(a1[1],a2[1])  
   
       def __stoppingcriterium(self,norm_r,r,p):  
           return self.stoppingcriterium(r[1],r[0],p)  
   
       def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
           return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)  
   
       def setTolerance(self,tolerance=1.e-8):  
               self.__tol=tolerance  
       def getTolerance(self):  
               return self.__tol  
       def setToleranceReductionFactor(self,reduction=0.01):  
               self.__reduction=reduction  
       def getSubProblemTolerance(self):  
               return self.__reduction*self.getTolerance()  
   
       def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):  
               """  
               solves the saddle point problem using initial guesses v and p.  
   
               @param max_iter: maximum number of iteration steps.  
               """  
               self.verbose=verbose  
               self.show_details=show_details and self.verbose  
   
               # assume p is known: then v=A^-1(f-B^*p)  
               # which leads to BA^-1B^*p = BA^-1f    
   
           # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)        
           self.__z=v+self.solve_A(v,p*0)  
               Bz=self.B(self.__z)  
               #  
           #   solve BA^-1B^*p = Bz  
               #  
               #  
               #  
               self.iter=0  
           if solver=='GMRES':        
                 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter  
                 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='TFQMR':        
                 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter  
                 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='MINRES':        
                 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter  
                 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
                 
           if solver=='GMRESC':        
                 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter  
                 p0=self.solve_prec1(Bz)  
             #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))  
                 #p-=alfa  
                 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)  
                 #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())  
   
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
             p=x[1]  
         u=v+self.solve_A(v,p)        
         #p=x[1]  
         #u=x[0]  
   
               if solver=='PCG':  
                 #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
                 #  
                 #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)  
                 #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)  
                 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter  
                 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)  
             u=r[0]    
                 # print "DIFF=",util.integrate(p)  
   
               # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)  
1502    
1503            return u,p        def __inner_PCG(self,p,r):
1504             return self.inner_pBv(p,r[1])
1505    
1506        def __Msolve(self,r):        def __Msolve_PCG(self,r):
1507            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1508          #=============================================================
1509          def __Aprod_GMRES(self,x):
1510              # return w,q  from v, p
1511              # solve Aw =Av+B^*p as Aw =f-A(z-v)-B^*(-p)
1512              #  and  Sq=B(v-w)
1513              v=x[0]
1514              p=x[1]
1515              w=self.solve_A(self.__z-v,-p)
1516              Bw=self.B(w-v)
1517          q=self.solve_prec(Bw)
1518              return ArithmeticTuple(w,q)
1519    
1520          def __inner_GMRES(self,a1,a2):
1521             return self.inner_p(a1[1],a2[1])+self.inner_v(a1[0],a2[0])
1522    
1523          #=============================================================
1524          def norm(self,v,p):
1525            f=self.inner_p(p,p)+self.inner_v(v,v)
1526            if f<0:
1527                raise ValueError,"negative norm."
1528            return math.sqrt(f)
1529    
1530        def __Msolve2(self,r):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):
1531            return self.solve_prec(r*1)           """
1532             solves the saddle point problem using initial guesses v and p.
1533    
1534        def __Mempty(self,r):           @param v: initial guess for velocity
1535            return r           @param p: initial guess for pressure
1536             @type v: L{Data}
1537             @type p: L{Data}
1538             @param useUzawa: indicate the usage of the Uzawa scheme. Otherwise the problem is solved in coupled form.
1539             @param max_iter: maximum number of iteration steps per correction attempt.
1540             @param verbose: show information on the progress of the saddlepoint problem solver.
1541             @param show_details: show details of the sub problems solves
1542             @param iter_restart: restart the iteration after C{iter_restart} steps (only used if useUzaw=False)
1543             @param max_correction_steps: maximum number of iteration steps in the attempt get |Bv| to zero.
1544             @return: new approximations for velocity and pressure
1545             @type useUzawa: C{bool}
1546             @type max_iter: C{int}
1547             @type verbose: C{bool}
1548             @type show_details: C{bool}
1549             @type iter_restart: C{int}
1550             @type max_correction_steps: C{int}
1551             @rtype: C{tuple} of L{Data} objects
1552             """
1553             self.verbose=verbose
1554             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)
1557             self.__z=v+self.solve_A(v,p*0)
1558             # tolerances:
1559             rtol=self.getTolerance()
1560             atol=self.getAbsoluteTolerance()
1561             if useUzawa:
1562                p0=self.solve_prec(self.B(self.__z))
1563                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
1568             if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1569             if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)
1570             # initialization
1571             self.iter=0
1572             correction_step=0
1573             converged=False
1574             # initial guess:
1575             q=p*1
1576             u=v*1
1577             while not converged :
1578                if useUzawa:
1579                   # assume p is known: then v=z-A^-1B^*p
1580                   #      
1581                   # which leads to BA^-1B^*p = Bz
1582                   #
1583                   # 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                   # 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())
1594                   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]  
1596               Bu=r[1]
1597                else:
1598                   #
1599                   #  with v=dv+z
1600                   #
1601                   #   A v + B p = f
1602                   #   B v       = 0
1603                   #
1604                   # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1})  -S^{-1}]]
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())
1608                   x=GMRES(ArithmeticTuple(w,self.solve_prec(self.B(u+w))),self.__Aprod_GMRES, ArithmeticTuple(u,q), \
1609                             self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1610               u=x[0]
1611                   q=x[1]
1612               Bu=self.B(u)
1613                # now we check if |Bu| ~ 0 or more precise |Bu|_p  <= rtol * |v|_v
1614                nu=self.inner_v(u,u)
1615                p2=self.solve_prec(Bu)
1616                nBu=self.inner_p(p2,p2)
1617                if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."
1618                nu=math.sqrt(nu)
1619                nBu=math.sqrt(nBu)
1620                if self.verbose: print "saddle point solver: norm v= %e (Bv = %e)"%(nu,nBu)
1621                QTOL=atol+nu*rtol
1622                if nBu <= QTOL:
1623                    converged=True
1624                else:
1625                    ATOL=QTOL/nBu*ATOL*0.3
1626                    if self.verbose: print "correction step %s: tolerance reduced to %e."%(correction_step,ATOL)
1627                    converged=False
1628                correction_step+=1
1629                if correction_step>max_correction_steps:
1630                   raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1631             if self.verbose: print "saddle point solver: tolerance reached."
1632         return u,q
1633    
1634          #==========================================================================================================================
1635          def setTolerance(self,tolerance=1.e-4):
1636             """
1637             sets the relative tolerance for (v,p)
1638    
1639        def __Aprod(self,p):           @param tolerance: tolerance to be used
1640            # return BA^-1B*p           @type tolerance: non-negative C{float}
1641            #solve Av =B^*p as Av =f-Az-B^*(-p)           """
1642            v=self.solve_A(self.__z,-p)           if tolerance<0:
1643            return ArithmeticTuple(v, self.B(v))               raise ValueError,"tolerance must be positive."
1644             self.__rtol=tolerance
1645             self.setSubProblemTolerance()
1646    
1647        def __Anext(self,x):        def getTolerance(self):
1648            # return next v,p           """
1649            #solve Av =-B^*p as Av =f-Az-B^*p           returns the relative tolerance
1650    
1651        pc=x[1]           @return: relative tolerance
1652            v=self.solve_A(self.__z,-pc)           @rtype: C{float}
1653        p=self.solve_prec1(self.B(v))           """
1654             return self.__rtol
1655          def setAbsoluteTolerance(self,tolerance=0.):
1656             """
1657             sets the absolute tolerance
1658    
1659            return ArithmeticTuple(v,p)           @param tolerance: tolerance to be used
1660             @type tolerance: non-negative C{float}
1661             """
1662             if tolerance<0:
1663                 raise ValueError,"tolerance must be non-negative."
1664             self.__atol=tolerance
1665          def getAbsoluteTolerance(self):
1666             """
1667             returns the absolute tolerance
1668    
1669             @return: absolute tolerance
1670             @rtype: C{float}
1671             """
1672             return self.__atol
1673    
1674        def __Aprod2(self,p):        def setSubProblemTolerance(self,rtol=None):
1675            # return BA^-1B*p           """
1676            #solve Av =B^*p as Av =f-Az-B^*(-p)           sets the relative tolerance to solve the subproblem(s).
       v=self.solve_A(self.__z,-p)  
           return self.B(v)  
1677    
1678        def __Aprod_Newton(self,p):           @param rtol: relative tolerence
1679            # return BA^-1B*p - Bz           @type rtol: positive C{float}
1680            #solve Av =-B^*p as Av =f-Az-B^*p           """
1681        v=self.solve_A(self.__z,-p)           if rtol == None:
1682            return self.B(v-self.__z)                rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1683             if rtol<=0:
1684        def __Aprod_Newton2(self,x):               raise ValueError,"tolerance must be positive."
1685            # return BA^-1B*p - Bz           self.__sub_tol=rtol
1686            #solve Av =-B^*p as Av =f-Az-B^*p        def getSubProblemTolerance(self):
1687            pc=x[1]           """
1688        v=self.solve_A(self.__z,-pc)           returns the subproblem reduction factor
           p=self.solve_prec1(self.B(v-self.__z))  
           return ArithmeticTuple(v,p)  
1689    
1690             @return: subproblem reduction factor
1691             @rtype: C{float}
1692             """
1693             return self.__sub_tol
1694    
1695  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1696     """     """
1697     creates a mask on the Solution(domain) function space which one for samples     create a mask on the Solution(domain) function space which one for samples
1698     that touch the boundary tagged by tags.     that touch the boundary tagged by tags.
1699    
1700     usage: m=MaskFromBoundaryTag(domain,"left", "right")     usage: m=MaskFromBoundaryTag(domain,"left", "right")

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

  ViewVC Help
Powered by ViewVC 1.1.26