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

revision 1469 by gross, Thu Apr 3 05:16:56 2008 UTC revision 1557 by artak, Mon May 19 04:44:27 2008 UTC
# Line 474  class IterationHistory(object): Line 474  class IterationHistory(object):
474
475         """         """
476         self.history.append(norm_r)         self.history.append(norm_r)
477         if self.verbose: print "iter: %s:  inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])         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]         return self.history[-1]<=self.tolerance * self.history[0]
479
480     def stoppingcriterium2(self,norm_r,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}         returns True if the C{norm_r} is C{tolerance}*C{norm_b}
483
# Line 490  class IterationHistory(object): Line 490  class IterationHistory(object):
490         @rtype: C{bool}         @rtype: C{bool}
491
492         """         """
493           if TOL==None:
494              TOL=self.tolerance
495         self.history.append(norm_r)         self.history.append(norm_r)
496         if self.verbose: print "iter: %s:  norm(r) = %e"%(len(self.history)-1, self.history[-1])         if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])
497         return self.history[-1]<=self.tolerance * norm_b         return self.history[-1]<=TOL * norm_b
498
499  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
500     """     """
# Line 579  def givapp(c,s,vin): Line 581  def givapp(c,s,vin):
581              vrot[i:i+2]=w1,w2              vrot[i:i+2]=w1,w2
582      return vrot      return vrot
583
584  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  ##############################################
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     iter=0
599     r=Msolve(b)     r=Msolve(b)
600     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
# Line 593  def GMRES(b, Aprod, Msolve, bilinearform Line 608  def GMRES(b, Aprod, Msolve, bilinearform
608        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
609        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm,"negative norm."
610
611     h=numarray.zeros((iter_max,iter_max),numarray.Float64)     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
612     c=numarray.zeros(iter_max,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
613     s=numarray.zeros(iter_max,numarray.Float64)     s=numarray.zeros(iter_restart,numarray.Float64)
614     g=numarray.zeros(iter_max,numarray.Float64)     g=numarray.zeros(iter_restart,numarray.Float64)
615     v=[]     v=[]
616
617     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
618
619     v.append(r/rho)     v.append(r/rho)
620     g[0]=rho     g[0]=rho
621
622     while not stoppingcriterium(rho,norm_b):     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
623
624      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
625

626      p=Msolve(Aprod(v[iter]))      p=Msolve(Aprod(v[iter]))
627
628      v.append(p)      v.append(p)
# Line 615  def GMRES(b, Aprod, Msolve, bilinearform Line 630  def GMRES(b, Aprod, Msolve, bilinearform
630      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
631
632  # Modified Gram-Schmidt  # Modified Gram-Schmidt
633        for j in range(iter+1):
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):      for j in range(iter+1):
774        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j][iter]=bilinearform(v[j],v[iter+1])
775        v[iter+1]+=(-1.)*h[j][iter]*v[j]        v[iter+1]+=(-1.)*h[j][iter]*v[j]
# Line 656  def GMRES(b, Aprod, Msolve, bilinearform Line 811  def GMRES(b, Aprod, Msolve, bilinearform
811          rho=abs(g[iter+1])          rho=abs(g[iter+1])
812      iter+=1      iter+=1
813
# At this point either iter > iter_max or rho < tol.
# It's time to compute x and leave.

814     if iter > 0 :     if iter > 0 :
815       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)
816       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1][iter-1]
# Line 673  def GMRES(b, Aprod, Msolve, bilinearform Line 825  def GMRES(b, Aprod, Msolve, bilinearform
825     else : xhat=v[0]     else : xhat=v[0]
826
827     x += xhat     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     return x          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
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):  class ArithmeticTuple(object):
# Line 729  class ArithmeticTuple(object): Line 1205  class ArithmeticTuple(object):
1205         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1206         """         """
1207         out=[]         out=[]
1208         for i in range(len(self)):         other=1.*other
1209           if isinstance(other,float):
1210        for i in range(len(self)):
1211             out.append(self[i]*other)             out.append(self[i]*other)
1212           else:
1213            for i in range(len(self)):
1214               out.append(self[i]*other[i])
1215         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1216
1217     def __rmul__(self,other):     def __rmul__(self,other):
# Line 743  class ArithmeticTuple(object): Line 1224  class ArithmeticTuple(object):
1224         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1225         """         """
1226         out=[]         out=[]
1227         for i in range(len(self)):         other=1.*other
1228           if isinstance(other,float):
1229        for i in range(len(self)):
1230             out.append(other*self[i])             out.append(other*self[i])
1231           else:
1232            for i in range(len(self)):
1233               out.append(other[i]*self[i])
1234         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1235
1236  #########################  #########################
# Line 760  class ArithmeticTuple(object): Line 1246  class ArithmeticTuple(object):
1246         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1247         """         """
1248         out=[]         out=[]
1249         for i in range(len(self)):         other=1.*other
1250             out.append(self[i]/other)         if isinstance(other,float):
1251        for i in range(len(self)):
1252               out.append(self[i]*(1./other))
1253           else:
1254            for i in range(len(self)):
1255               out.append(self[i]*(1./other[i]))
1256         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1257
1258     def __rdiv__(self,other):     def __rdiv__(self,other):
# Line 774  class ArithmeticTuple(object): Line 1265  class ArithmeticTuple(object):
1265         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1266         """         """
1267         out=[]         out=[]
1268         for i in range(len(self)):         other=1.*other
1269             out.append(other/self[i])         if isinstance(other,float):
1270            for i in range(len(self)):
1271               out.append(other*(1./self[i]))
1272           else:
1273            for i in range(len(self)):
1274               out.append(other[i]*(1./self[i]))
1275         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1276
1277  ##########################################33  ##########################################33
# Line 793  class ArithmeticTuple(object): Line 1289  class ArithmeticTuple(object):
1289             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1290         return self         return self
1291
1293           """
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
1347        """        """
1348        This provides a framwork for solving homogeneous saddle point problem of the form        This provides a framwork for solving homogeneous saddle point problem of the form
1414
1415        def __inner_p(self,p1,p2):        def __inner_p(self,p1,p2):
1416           return self.inner(p1,p2)           return self.inner(p1,p2)
1417
1418          def __inner_a(self,a1,a2):
1419             return self.inner_a(a1,a2)
1420
1421          def __inner_a1(self,a1,a2):
1422             return self.inner(a1[1],a2[1])
1423
1424        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
1425            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
1426
1427        def __stoppingcriterium_GMRES(self,norm_r,norm_b):        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1428            return self.stoppingcriterium_GMRES(norm_r,norm_b)            return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1429
1430        def setTolerance(self,tolerance=1.e-8):        def setTolerance(self,tolerance=1.e-8):
1431                self.__tol=tolerance                self.__tol=tolerance
1436        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1437                return self.__reduction*self.getTolerance()                return self.__reduction*self.getTolerance()
1438
1439        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):        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.                solves the saddle point problem using initial guesses v and p.
1442
1449                # which leads to BA^-1B^*p = BA^-1f                  # which leads to BA^-1B^*p = BA^-1f
1450
1451            # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)                  # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)

1452            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
1453                  Bz=self.B(self.__z)
Bz=self.B(self.__z)
1454                #                #
1455            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1456                #                #
#   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1457                #                #
#   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
#                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1458                #                #
1459                self.iter=0                self.iter=0
1460            if solver=='GMRES':                  if solver=='GMRES':
1461                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1462                  p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)                  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                  # solve Au=f-B^*p
1464                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1465                  #       u=v+(u-v)                  #       u=v+(u-v)
1466          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1467
1468                else:            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                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1515                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1516              u=r[0]                u=r[0]
1517                print "div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                  print "DIFF=",util.integrate(p)
1518
1519                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1520
1521            return u,p            return u,p
1522
1523        def __Msolve(self,r):        def __Msolve(self,r):
1524            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1525
1526        def __Msolve_GMRES(self,r):        def __Msolve2(self,r):
1527            return self.solve_prec(r)            return self.solve_prec(r*1)
1528
1529          def __Mempty(self,r):
1530              return r
1531
1532
1533        def __Aprod(self,p):        def __Aprod(self,p):
1534            # return BA^-1B*p            # return BA^-1B*p
1535            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1536            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1537            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1538
1539        def __Aprod_GMRES(self,p):        def __Anext(self,x):
1540            # return BA^-1B*p            # return next v,p
1541            #solve Av =-B^*p as Av =f-Az-B^*p            #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)        v=self.solve_A(self.__z,-p)
1554            return self.B(v)            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
1571     """     """
1572     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem