# Diff of /trunk/escript/py_src/pdetools.py

revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1557 by artak, Mon May 19 04:44:27 2008 UTC
# Line 48  import numarray Line 48  import numarray
48  import util  import util
49  import math  import math
50
51    ##### Added by Artak
52    # from Numeric import zeros,Int,Float64
53    ###################################
54
55
56  class TimeIntegrationManager:  class TimeIntegrationManager:
57    """    """
58    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
# Line 435  class NegativeNorm(SolverSchemeException Line 440  class NegativeNorm(SolverSchemeException
440     """     """
441     pass     pass
442
443  def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)):  class IterationHistory(object):
444     """     """
445     Solver for     The IterationHistory class is used to define a stopping criterium. It keeps track of the
446       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
447       a given tolerance.
448       """
449       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
450          """
451          Initialization
452
453     M{Ax=b}        @param tolerance: tolerance
454          @type tolerance: positive C{float}
455          @param verbose: switches on the printing out some information
456          @type verbose: C{bool}
457          """
458          if not tolerance>0.:
459              raise ValueError,"tolerance needs to be positive."
460          self.tolerance=tolerance
461          self.verbose=verbose
462          self.history=[]
463       def stoppingcriterium(self,norm_r,r,x):
464           """
465           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.
466
467     with a symmetric and positive definite operator A (more details required!).
468     It uses the conjugate gradient method with preconditioner M providing an approximation of A.         @param norm_r: current residual norm
469           @type norm_r: non-negative C{float}
470           @param r: current residual (not used)
471           @param x: current solution approximation (not used)
472           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
473           @rtype: C{bool}
474
475     The iteration is terminated if         """
476           self.history.append(norm_r)
477           if self.verbose: print "iter: #s:  inner(rhat,r) = #e"#(len(self.history)-1, self.history[-1])
478           return self.history[-1]<=self.tolerance * self.history[0]
479
480     M{norm(r) <= tolerance * norm(b)}     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
481           """
482           returns True if the C{norm_r} is C{tolerance}*C{norm_b}
483
484     where C{norm()} defines a norm and
485           @param norm_r: current residual norm
486           @type norm_r: non-negative C{float}
487           @param norm_b: norm of right hand side
488           @type norm_b: non-negative C{float}
489           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
490           @rtype: C{bool}
491
492     M{r = b-Ax}         """
493           if TOL==None:
494              TOL=self.tolerance
495           self.history.append(norm_r)
496           if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])
497           return self.history[-1]<=TOL * norm_b
498
499    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
500       """
501       Solver for
502
503     the residual.     M{Ax=b}
504
505       with a symmetric and positive definite operator A (more details required!).
506       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
507
508       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
509
510     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
511
# Line 461  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 514  def PCG(b,x,Aprod,Msolve,bilinearform, n
514     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst.
515
516     @param b: the right hand side of the liner system. C{b} is altered.     @param b: the right hand side of the liner system. C{b} is altered.
517     @type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y)     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
@param x: an initial guess for the solution
@type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y)
518     @param Aprod: returns the value Ax     @param Aprod: returns the value Ax
519     @type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R.     @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}.
520     @param Msolve: solves Mx=r     @param Msolve: solves Mx=r
521     @type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S.     @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
522    type like argument C{x}.
523     @param bilinearform: inner product C{<x,r>}     @param bilinearform: inner product C{<x,r>}
524     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}.     @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}.
525     @param norm: norm calculation for the residual C{r=b-Ax}.     @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.
526     @type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}.     @type stoppingcriterium: function that returns C{True} or C{False}
527     @param verbose: switches on the printing out some information     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
528     @type verbose: C{bool}     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
529     @param iter_max: maximum number of iteration steps.     @param iter_max: maximum number of iteration steps.
530     @type iter_max: C{int}     @type iter_max: C{int}
531     @param tolerance: tolerance     @return: the solution approximation and the corresponding residual
532     @type tolerance: positive C{float}     @rtype: C{tuple}
533     @return: the solution apprximation and the corresponding residual     @warning: C{b} and C{x} are altered.
@rtype: C{tuple} of an S type and and an R type object.A
@warning: C{b} ans C{x} are altered.
534     """     """
if verbose:
print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance   =%e"%(iter_max,tolerance)
535     iter=0     iter=0
536     normb = norm(b)     if x==None:
537     if normb<0: raise NegativeNorm        x=0*b
538       else:
539     b += (-1)*Aprod(x)        b += (-1)*Aprod(x)
540     r=b     r=b
541     rhat=Msolve(r)     rhat=Msolve(r)
542     d = rhat;     d = rhat
543     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
544       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
545
546     while True:     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
547         normr=norm(r)         iter+=1
if normr<0: raise NegativeNorm
if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb)
if normr <= tolerance * normb: return x,r

iter+=1 # k = iter = 1 first time through
548         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
549
550         q=Aprod(d)         q=Aprod(d)
# Line 515  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 559  def PCG(b,x,Aprod,Msolve,bilinearform, n
559         d=rhat         d=rhat
560
561         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
562           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
563
564       return x,r
565
566
567    ############################
568    # Added by Artak
569    #################################3
570
571    #Apply a sequence of k Givens rotations, used within gmres codes
572    # vrot=givapp(c, s, vin, k)
573    def givapp(c,s,vin):
574        vrot=vin # warning: vin is altered!!!!
575        if isinstance(c,float):
576            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
577        else:
578            for i in range(len(c)):
579                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
580            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
581                vrot[i:i+2]=w1,w2
582        return vrot
583
584    ##############################################
585    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
586    ################################################
587       m=iter_restart
588       iter=0
589       xc=x
590       while True:
591          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
592          xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
593          if stopped: break
594          iter+=iter_restart
595       return xc
596
597    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
598       iter=0
599       r=Msolve(b)
600       r_dot_r = bilinearform(r, r)
601       if r_dot_r<0: raise NegativeNorm,"negative norm."
602       norm_b=math.sqrt(r_dot_r)
603
604       if x==None:
605          x=0*b
606       else:
607          r=Msolve(b-Aprod(x))
608          r_dot_r = bilinearform(r, r)
609          if r_dot_r<0: raise NegativeNorm,"negative norm."
610
611       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
612       c=numarray.zeros(iter_restart,numarray.Float64)
613       s=numarray.zeros(iter_restart,numarray.Float64)
614       g=numarray.zeros(iter_restart,numarray.Float64)
615       v=[]
616
617       rho=math.sqrt(r_dot_r)
618
619       v.append(r/rho)
620       g[0]=rho
621
622       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
625
626        p=Msolve(Aprod(v[iter]))
627
628        v.append(p)
629
630        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
631
632    # Modified Gram-Schmidt
633        for j in range(iter+1):
634          h[j][iter]=bilinearform(v[j],v[iter+1])
635          v[iter+1]-=h[j][iter]*v[j]
636
637        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
638        v_norm2=h[iter+1][iter]
639
640    # Reorthogonalize if needed
641        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
642         for j in range(iter+1):
643            hr=bilinearform(v[j],v[iter+1])
644                h[j][iter]=h[j][iter]+hr
645                v[iter+1] -= hr*v[j]
646
647         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
648         h[iter+1][iter]=v_norm2
649
650    #   watch out for happy breakdown
651            if not v_norm2 == 0:
652             v[iter+1]=v[iter+1]/h[iter+1][iter]
653
654    #   Form and store the information for the new Givens rotation
655        if iter > 0 :
656            hhat=numarray.zeros(iter+1,numarray.Float64)
657            for i in range(iter+1) : hhat[i]=h[i][iter]
658            hhat=givapp(c[0:iter],s[0:iter],hhat);
659                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])
662
663        if mu!=0 :
664            c[iter]=h[iter][iter]/mu
665            s[iter]=-h[iter+1][iter]/mu
666            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
667            h[iter+1][iter]=0.0
668            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
669
670    # Update the residual norm
671
672            rho=abs(g[iter+1])
673        iter+=1
674
675    # At this point either iter > iter_max or rho < tol.
676    # It's time to compute x and leave.
677
678       if iter > 0 :
679         y=numarray.zeros(iter,numarray.Float64)
680         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
681         if iter > 1 :
682            i=iter-2
683            while i>=0 :
684              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
685              i=i-1
686         xhat=v[iter-1]*y[iter-1]
687         for i in range(iter-1):
688        xhat += v[i]*y[i]
689       else : xhat=v[0]
690
691       x += xhat
692       if iter<iter_restart-1:
693          stopped=True
694       else:
695          stopped=False
696
697       return x,stopped
698
699
700    ######################################################
701    def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
702    ######################################################
703
704    # DIRDER estimates the directional derivative of a function F.
705
706
707    # Hardwired difference increment.
708    #
709      epsnew = 1.0e-07
710    #
711    #  Scale the step.
712    #
713      norm_w=math.sqrt(bilinearform(w,w))
714      if norm_w== 0.0:
715        return x/x
716
717      epsnew = epsnew / norm_w
718
719      if norm_w > 0.0:
720        epsnew = epsnew * math.sqrt(bilinearform(x,x))
721    #
722    #  DEL and F1 could share the same space if storage
723    #  is more important than clarity.
724    #
725
726      DEL = x + epsnew * w
727      f1 = -Msolve(Aprod(DEL))
728      z = ( f1 - f0 ) / epsnew
729      return z
730
731    ######################################################
732    def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
733    ######################################################
734       b=-f0
735       b_dot_b = bilinearform(b, b)
736       if b_dot_b<0: raise NegativeNorm,"negative norm."
737       norm_b=math.sqrt(b_dot_b)
738
739       r=b
740
741       if x==None:
742          x=0*f0
743       else:
744          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0
745
746       r_dot_r = bilinearform(r, r)
747       if r_dot_r<0: raise NegativeNorm,"negative norm."
748
749       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
750       c=numarray.zeros(iter_restart,numarray.Float64)
751       s=numarray.zeros(iter_restart,numarray.Float64)
752       g=numarray.zeros(iter_restart,numarray.Float64)
753       v=[]
754
755       rho=math.sqrt(r_dot_r)
756
757       v.append(r/rho)
758       g[0]=rho
759       iter=0
760
761       while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
762
763        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
764
765
766            p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
767
768        v.append(p)
769
770        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
771
772    # Modified Gram-Schmidt
773        for j in range(iter+1):
774          h[j][iter]=bilinearform(v[j],v[iter+1])
775          v[iter+1]+=(-1.)*h[j][iter]*v[j]
776
777        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
778        v_norm2=h[iter+1][iter]
779
780
781    # Reorthogonalize if needed
782        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
783         for j in range(iter+1):
784            hr=bilinearform(v[j],v[iter+1])
785                h[j][iter]=h[j][iter]+hr #vhat
786                v[iter+1] +=(-1.)*hr*v[j]
787
788         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
789         h[iter+1][iter]=v_norm2
790
791    #   watch out for happy breakdown
792            if v_norm2 != 0:
793             v[iter+1]=v[iter+1]/h[iter+1][iter]
794
795    #   Form and store the information for the new Givens rotation
796        if iter > 0 :
797            hhat=[]
798            for i in range(iter+1) : hhat.append(h[i][iter])
799            hhat=givapp(c[0:iter],s[0:iter],hhat);
800                for i in range(iter+1) : h[i][iter]=hhat[i]
801
802        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
803        if mu!=0 :
804            c[iter]=h[iter][iter]/mu
805            s[iter]=-h[iter+1][iter]/mu
806            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
807            h[iter+1][iter]=0.0
808            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
809
810    # Update the residual norm
811            rho=abs(g[iter+1])
812        iter+=1
813
814       if iter > 0 :
815         y=numarray.zeros(iter,numarray.Float64)
816         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
817         if iter > 1 :
818            i=iter-2
819            while i>=0 :
820              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
821              i=i-1
822         xhat=v[iter-1]*y[iter-1]
823         for i in range(iter-1):
824        xhat += v[i]*y[i]
825       else : xhat=v[0]
826
827       x += xhat
828       if iter<iter_restart-1:
829          stopped=True
830       else:
831          stopped=False
832
833       return x,stopped
834
835    #################################################
836    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
837    #################################################
838        #
839        #  minres solves the system of linear equations Ax = b
840        #  where A is a symmetric matrix (possibly indefinite or singular)
841        #  and b is a given vector.
842        #
843        #  "A" may be a dense or sparse matrix (preferably sparse!)
844        #  or the name of a function such that
845        #               y = A(x)
846        #  returns the product y = Ax for any given vector x.
847        #
848        #  "M" defines a positive-definite preconditioner M = C C'.
849        #  "M" may be a dense or sparse matrix (preferably sparse!)
850        #  or the name of a function such that
851        #  solves the system My = x for any given vector x.
852        #
853        #
854
855        #------------------------------------------------------------------
856        # Set up y and v for the first Lanczos vector v1.
857        # y  =  beta1 P' v1,  where  P = C**(-1).
858        # v is really P' v1.
859        #------------------------------------------------------------------
860        if x==None:
861          x=0*b
862        else:
863          b += (-1)*Aprod(x)
864
865        r1    = b
866        y = Msolve(b)
867        beta1 = bilinearform(y,b)
868
869        if beta1< 0: raise NegativeNorm,"negative norm."
870
871        #  If b = 0 exactly, stop with x = 0.
872        if beta1==0: return x*0.
873
874        if beta1> 0:
875          beta1  = math.sqrt(beta1)
876
877        #------------------------------------------------------------------
878        # Initialize quantities.
879        # ------------------------------------------------------------------
880        iter   = 0
881        Anorm = 0
882        ynorm = 0
883        oldb   = 0
884        beta   = beta1
885        dbar   = 0
886        epsln  = 0
887        phibar = beta1
888        rhs1   = beta1
889        rhs2   = 0
890        rnorm  = phibar
891        tnorm2 = 0
892        ynorm2 = 0
893        cs     = -1
894        sn     = 0
895        w      = b*0.
896        w2     = b*0.
897        r2     = r1
898        eps    = 0.0001
899
900        #---------------------------------------------------------------------
901        # Main iteration loop.
902        # --------------------------------------------------------------------
903        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
904
905        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
906            iter    = iter  +  1
907
908            #-----------------------------------------------------------------
909            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
910            # The general iteration is similar to the case k = 1 with v0 = 0:
911            #
912            #   p1      = Operator * v1  -  beta1 * v0,
913            #   alpha1  = v1'p1,
914            #   q2      = p2  -  alpha1 * v1,
915            #   beta2^2 = q2'q2,
916            #   v2      = (1/beta2) q2.
917            #
918            # Again, y = betak P vk,  where  P = C**(-1).
919            #-----------------------------------------------------------------
920            s = 1/beta                 # Normalize previous vector (in y).
921            v = s*y                    # v = vk if P = I
922
923            y      = Aprod(v)
924
925            if iter >= 2:
926              y = y - (beta/oldb)*r1
927
928            alfa   = bilinearform(v,y)              # alphak
929            y      += (- alfa/beta)*r2
930            r1     = r2
931            r2     = y
932            y = Msolve(r2)
933            oldb   = beta                           # oldb = betak
934            beta   = bilinearform(y,r2)             # beta = betak+1^2
935            if beta < 0: raise NegativeNorm,"negative norm."
936
937            beta   = math.sqrt( beta )
938            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
939
940            if iter==1:                 # Initialize a few things.
941              gmax   = abs( alfa )      # alpha1
942              gmin   = gmax             # alpha1
943
944            # Apply previous rotation Qk-1 to get
945            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
946            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
947
948            oldeps = epsln
949            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
950            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
951            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
952            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
953
954            # Compute the next plane rotation Qk
955
956            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
957            gamma  = max(gamma,eps)
958            cs     = gbar / gamma             # ck
959            sn     = beta / gamma             # sk
960            phi    = cs * phibar              # phik
961            phibar = sn * phibar              # phibark+1
962
963            # Update  x.
964
965            denom = 1/gamma
966            w1    = w2
967            w2    = w
968            w     = (v - oldeps*w1 - delta*w2) * denom
969            x     +=  phi*w
970
971            # Go round again.
972
973            gmax   = max(gmax,gamma)
974            gmin   = min(gmin,gamma)
975            z      = rhs1 / gamma
976            ynorm2 = z*z  +  ynorm2
977            rhs1   = rhs2 -  delta*z
978            rhs2   =      -  epsln*z
979
980            # Estimate various norms and test for convergence.
981
982            Anorm  = math.sqrt( tnorm2 )
983            ynorm  = math.sqrt( ynorm2 )
984
985            rnorm  = phibar
986
987        return x
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
993        lmaxit=100
994        etamax=.5
995
996        n = 1 #len(x)
997        iter=0
998
999        # evaluate f at the initial iterate
1000        # compute the stop tolerance
1001        #
1002        r=b
1003        if x==None:
1004          x=0*b
1005        else:
1006          r += (-1)*Aprod(x)
1007
1008        f0=-Msolve(r)
1009        fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1010        fnrmo=1
1011        stop_tol=atol + rtol*fnrm
1012        #
1013        # main iteration loop
1014        #
1015        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
1018            #
1019            # keep track of the ratio (rat = fnrm/frnmo)
1020            # of successive residual norms and
1021            # the iteration counter (iter)
1022            #
1023            #rat=fnrm/fnrmo
1024            fnrmo=fnrm
1025            iter+=1
1026            #
1027                # compute the step using a GMRES(m) routine especially designed for this purpose
1028            #
1029                initer=0
1030                while True:
1031                   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
1033                   initer+=iter_restart
1034            x+=xc
1035            f0=-Msolve(Aprod(x))
1036            fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1037            rat=fnrm/fnrmo
1038
1039
1040        #   adjust eta
1041        #
1042            if etamax > 0:
1043                etaold=etamax
1044                etanew=gamma*rat*rat
1045                if gamma*etaold*etaold > .1 :
1046                    etanew=max(etanew,gamma*etaold*etaold)
1047
1048                etamax=min(etanew,etamax)
1049                etamax=max(etamax,.5*stop_tol/fnrm)
1050
1051        return x
1052
1053    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1054
1055    # TFQMR solver for linear systems
1056    #
1057    #
1058    # initialization
1059    #
1060      errtol = math.sqrt(bilinearform(b,b))
1061      norm_b=errtol
1062      kmax  = iter_max
1063      error = []
1064
1065      if math.sqrt(bilinearform(x,x)) != 0.0:
1066        r = b - Aprod(x)
1067      else:
1068        r = b
1069
1070      r=Msolve(r)
1071
1072      u1=0
1073      u2=0
1074      y1=0
1075      y2=0
1076
1077      w = r
1078      y1 = r
1079      iter = 0
1080      d = 0
1081
1082      v = Msolve(Aprod(y1))
1083      u1 = v
1084
1085      theta = 0.0;
1086      eta = 0.0;
1087      tau = math.sqrt(bilinearform(r,r))
1088      error = [ error, tau ]
1089      rho = tau * tau
1090      m=1
1091    #
1092    #  TFQMR iteration
1093    #
1094    #  while ( iter < kmax-1 ):
1095
1096      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1097        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1098
1099        sigma = bilinearform(r,v)
1100
1101        if ( sigma == 0.0 ):
1102          raise 'TFQMR breakdown, sigma=0'
1103
1104
1105        alpha = rho / sigma
1106
1107        for j in range(2):
1108    #
1109    #   Compute y2 and u2 only if you have to
1110    #
1111          if ( j == 1 ):
1112            y2 = y1 - alpha * v
1113            u2 = Msolve(Aprod(y2))
1114
1115          m = 2 * (iter+1) - 2 + (j+1)
1116          if j==0:
1117             w = w - alpha * u1
1118             d = y1 + ( theta * theta * eta / alpha ) * d
1119          if j==1:
1120             w = w - alpha * u2
1121             d = y2 + ( theta * theta * eta / alpha ) * d
1122
1123          theta = math.sqrt(bilinearform(w,w))/ tau
1124          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1125          tau = tau * theta * c
1126          eta = c * c * alpha
1127          x = x + eta * d
1128    #
1129    #  Try to terminate the iteration at each pass through the loop
1130    #
1131         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1132         #   error = [ error, tau ]
1133         #   total_iters = iter
1134         #   break
1135
1136
1137        if ( rho == 0.0 ):
1138          raise 'TFQMR breakdown, rho=0'
1139
1140
1141        rhon = bilinearform(r,w)
1142        beta = rhon / rho;
1143        rho = rhon;
1144        y1 = w + beta * y2;
1145        u1 = Msolve(Aprod(y1))
1146        v = u1 + beta * ( u2 + beta * v )
1147        error = [ error, tau ]
1148        total_iters = iter
1149
1150        iter = iter + 1
1151
1152      return x
1153
1154
1155    #############################################
1156
1157    class ArithmeticTuple(object):
1158       """
1159       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1160
1161       example of usage:
1162
1163       from esys.escript import Data
1164       from numarray import array
1165       a=Data(...)
1166       b=array([1.,4.])
1167       x=ArithmeticTuple(a,b)
1168       y=5.*x
1169
1170       """
1171       def __init__(self,*args):
1172           """
1173           initialize object with elements args.
1174
1175           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1176           """
1177           self.__items=list(args)
1178
1179       def __len__(self):
1180           """
1181           number of items
1182
1183           @return: number of items
1184           @rtype: C{int}
1185           """
1186           return len(self.__items)
1187
1188       def __getitem__(self,index):
1189           """
1190           get an item
1191
1192           @param index: item to be returned
1193           @type index: C{int}
1194           @return: item with index C{index}
1195           """
1196           return self.__items.__getitem__(index)
1197
1198       def __mul__(self,other):
1199           """
1200           scaling from the right
1201
1202           @param other: scaling factor
1203           @type other: C{float}
1204           @return: itemwise self*other
1205           @rtype: L{ArithmeticTuple}
1206           """
1207           out=[]
1208           other=1.*other
1209           if isinstance(other,float):
1210        for i in range(len(self)):
1211               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))
1216
1217       def __rmul__(self,other):
1218           """
1219           scaling from the left
1220
1221           @param other: scaling factor
1222           @type other: C{float}
1223           @return: itemwise other*self
1224           @rtype: L{ArithmeticTuple}
1225           """
1226           out=[]
1227           other=1.*other
1228           if isinstance(other,float):
1229        for i in range(len(self)):
1230               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))
1235
1236    #########################
1237    # Added by Artak
1238    #########################
1239       def __div__(self,other):
1240           """
1241           dividing from the right
1242
1243           @param other: scaling factor
1244           @type other: C{float}
1245           @return: itemwise self/other
1246           @rtype: L{ArithmeticTuple}
1247           """
1248           out=[]
1249           other=1.*other
1250           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))
1257
1258       def __rdiv__(self,other):
1259           """
1260           dividing from the left
1261
1262           @param other: scaling factor
1263           @type other: C{float}
1264           @return: itemwise other/self
1265           @rtype: L{ArithmeticTuple}
1266           """
1267           out=[]
1268           other=1.*other
1269           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))
1276
1277    ##########################################33
1278
1279       def __iadd__(self,other):
1280           """
1281           in-place add of other to self
1282
1283           @param other: increment
1284           @type other: C{ArithmeticTuple}
1285           """
1286           if len(self) != len(other):
1287               raise ValueError,"tuple length must match."
1288           for i in range(len(self)):
1289               self.__items[i]+=other[i]
1290           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):
1347          """
1348          This provides a framwork for solving homogeneous saddle point problem of the form
1349
1350                 Av+B^*p=f
1351                 Bv    =0
1352
1353          for the unknowns v and p and given operators A and B and given right hand side f.
1354          B^* is the adjoint operator of B is the given inner product.
1355
1356          """
1357          def __init__(self,**kwargs):
1358            self.setTolerance()
1359            self.setToleranceReductionFactor()
1360
1361          def initialize(self):
1362            """
1363            initialize the problem (overwrite)
1364            """
1365            pass
1366          def B(self,v):
1367             """
1368             returns Bv (overwrite)
1369             @rtype: equal to the type of p
1370
1371             @note: boundary conditions on p should be zero!
1372             """
1373             pass
1374
1375          def inner(self,p0,p1):
1376             """
1377             returns inner product of two element p0 and p1  (overwrite)
1378
1379             @type p0: equal to the type of p
1380             @type p1: equal to the type of p
1381             @rtype: C{float}
1382
1383             @rtype: equal to the type of p
1384             """
1385             pass
1386
1387          def solve_A(self,u,p):
1388             """
1389             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1390
1391             @rtype: equal to the type of v
1392             @note: boundary conditions on v should be zero!
1393             """
1394             pass
1395
1396          def solve_prec(self,p):
1397             """
1398             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1399
1400             @rtype: equal to the type of p
1401             """
1402             pass
1403
1404          def stoppingcriterium(self,Bv,v,p):
1405             """
1406             returns a True if iteration is terminated. (overwrite)
1407
1408             @rtype: C{bool}
1409             """
1410             pass
1411
1412          def __inner(self,p,r):
1413             return self.inner(p,r[1])
1414
1415          def __inner_p(self,p1,p2):
1416             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):
1425              return self.stoppingcriterium(r[1],r[0],p)
1426
1427          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1428              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1429
1430          def setTolerance(self,tolerance=1.e-8):
1431                  self.__tol=tolerance
1432          def getTolerance(self):
1433                  return self.__tol
1434          def setToleranceReductionFactor(self,reduction=0.01):
1435                  self.__reduction=reduction
1436          def getSubProblemTolerance(self):
1437                  return self.__reduction*self.getTolerance()
1438
1439          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1440                  """
1441                  solves the saddle point problem using initial guesses v and p.
1442
1443                  @param max_iter: maximum number of iteration steps.
1444                  """
1445                  self.verbose=verbose
1446                  self.show_details=show_details and self.verbose
1447
1448                  # assume p is known: then v=A^-1(f-B^*p)
1449                  # 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)
1452              self.__z=v+self.solve_A(v,p*0)
1453                  Bz=self.B(self.__z)
1454                  #
1455              #   solve BA^-1B^*p = Bz
1456                  #
1457                  #
1458                  #
1459                  self.iter=0
1460              if solver=='GMRES':
1461                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1462                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1463                    # solve Au=f-B^*p
1464                    #       A(u-v)=f-B^*p-Av
1465                    #       u=v+(u-v)
1466            u=v+self.solve_A(v,p)
1467
1468              if solver=='NewtonGMRES':
1469                    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.,atol=0,rtol=self.getTolerance())
1471                    # solve Au=f-B^*p
1472                    #       A(u-v)=f-B^*p-Av
1473                    #       u=v+(u-v)
1474            u=v+self.solve_A(v,p)
1475
1476
1477              if solver=='TFQMR':
1478                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1479                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1480                    # solve Au=f-B^*p
1481                    #       A(u-v)=f-B^*p-Av
1482                    #       u=v+(u-v)
1483            u=v+self.solve_A(v,p)
1484
1485              if solver=='MINRES':
1486                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1487                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1488                    # solve Au=f-B^*p
1489                    #       A(u-v)=f-B^*p-Av
1490                    #       u=v+(u-v)
1491            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':
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
1515                    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]
1517                    print "DIFF=",util.integrate(p)
1518
1519                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1520
1521              return u,p
1522
1523          def __Msolve(self,r):
1524              return self.solve_prec(r[1])
1525
1526          def __Msolve2(self,r):
1527              return self.solve_prec(r*1)
1528
1529          def __Mempty(self,r):
1530              return r
1531
1532
1533          def __Aprod(self,p):
1534              # return BA^-1B*p
1535              #solve Av =B^*p as Av =f-Az-B^*(-p)
1536              v=self.solve_A(self.__z,-p)
1537              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):
1551              # return BA^-1B*p
1552              #solve Av =B^*p as Av =f-Az-B^*(-p)
1553          v=self.solve_A(self.__z,-p)
1554              return self.B(v)
1555
1556          def __Aprod_Newton(self,p):
1557              # return BA^-1B*p - Bz
1558              #solve Av =-B^*p as Av =f-Az-B^*p
1559          v=self.solve_A(self.__z,-p)
1560              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     """     """
# Line 554  class SaddlePointProblem(object): Line 1604  class SaddlePointProblem(object):
1604         @param text: a text message         @param text: a text message
1605         @type text: C{str}         @type text: C{str}
1606         """         """
1607         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "#s: #s"%(str(self),text)
1608
1609     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1610         """         """
# Line 771  def MaskFromBoundaryTag(function_space,* Line 1821  def MaskFromBoundaryTag(function_space,*
1821     else:     else:
1822        return util.whereNonZero(util.interpolate(out,function_space))        return util.whereNonZero(util.interpolate(out,function_space))
1823
1824
1825

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

 ViewVC Help Powered by ViewVC 1.1.26