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

revision 1486 by artak, Thu Apr 10 00:10:50 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, iter_restart=10):  ##############################################
585    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
586    ################################################
587     m=iter_restart     m=iter_restart
588     iter=0     iter=0
589       xc=x
590     while True:     while True:
591        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
592        x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)        xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
iter+=iter_restart
593        if stopped: break        if stopped: break
594     return x        iter+=iter_restart
595       return xc
596
597  def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):  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 610  def GMRESm(b, Aprod, Msolve, bilinearfor Line 615  def GMRESm(b, Aprod, Msolve, bilinearfor
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
# Line 617  def GMRESm(b, Aprod, Msolve, bilinearfor Line 623  def GMRESm(b, Aprod, Msolve, bilinearfor
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 625  def GMRESm(b, Aprod, Msolve, bilinearfor Line 630  def GMRESm(b, Aprod, Msolve, bilinearfor
630      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
631
632  # Modified Gram-Schmidt  # Modified Gram-Schmidt
633        for j in range(iter+1):
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 666  def GMRESm(b, Aprod, Msolve, bilinearfor Line 811  def GMRESm(b, Aprod, Msolve, bilinearfor
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 683  def GMRESm(b, Aprod, Msolve, bilinearfor Line 825  def GMRESm(b, Aprod, Msolve, bilinearfor
825     else : xhat=v[0]     else : xhat=v[0]
826
827     x += xhat     x += xhat
828     if iter!=iter_restart-1:     if iter<iter_restart-1:
829        stopped=True        stopped=True
830     else:     else:
831        stopped=False        stopped=False
832
833     return x,stopped     return x,stopped
834
835    #################################################
836  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
837    #################################################
838      #      #
839      #  minres solves the system of linear equations Ax = b      #  minres solves the system of linear equations Ax = b
840      #  where A is a symmetric matrix (possibly indefinite or singular)      #  where A is a symmetric matrix (possibly indefinite or singular)
# Line 721  def MINRES(b, Aprod, Msolve, bilinearfor Line 864  def MINRES(b, Aprod, Msolve, bilinearfor
864
865      r1    = b      r1    = b
866      y = Msolve(b)      y = Msolve(b)
867      beta1 = bilinearform(b,y)      beta1 = bilinearform(y,b)
868
869      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
870
# Line 757  def MINRES(b, Aprod, Msolve, bilinearfor Line 900  def MINRES(b, Aprod, Msolve, bilinearfor
900      #---------------------------------------------------------------------      #---------------------------------------------------------------------
901      # Main iteration loop.      # Main iteration loop.
902      # --------------------------------------------------------------------      # --------------------------------------------------------------------
903      while not stoppingcriterium(rnorm,Anorm*ynorm):    #  checks ||r|| < (||A|| ||x||) * TOL      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      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
906          iter    = iter  +  1          iter    = iter  +  1
# Line 783  def MINRES(b, Aprod, Msolve, bilinearfor Line 926  def MINRES(b, Aprod, Msolve, bilinearfor
926            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
927
928          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
929          y      = (- alfa/beta)*r2 + y          y      += (- alfa/beta)*r2
930          r1     = r2          r1     = r2
931          r2     = y          r2     = y
932          y = Msolve(r2)          y = Msolve(r2)
933          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
934          beta   = bilinearform(r2,y)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
935          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm,"negative norm."
936
937          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
# Line 823  def MINRES(b, Aprod, Msolve, bilinearfor Line 966  def MINRES(b, Aprod, Msolve, bilinearfor
966          w1    = w2          w1    = w2
967          w2    = w          w2    = w
968          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
969          x     = x  +  phi*w          x     +=  phi*w
970
971          # Go round again.          # Go round again.
972
# Line 842  def MINRES(b, Aprod, Msolve, bilinearfor Line 985  def MINRES(b, Aprod, Msolve, bilinearfor
985          rnorm  = phibar          rnorm  = phibar
986
987      return x      return x
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 896  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 910  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 927  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 941  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 960  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)

def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
return self.stoppingcriterium_MINRES(norm_r,norm_Ax)

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='PCG'):        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
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                  # solve Au=f-B^*p
1481                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1482                  #       u=v+(u-v)                  #       u=v+(u-v)
1484
1485            if solver=='MINRES':                  if solver=='MINRES':
1486                  if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1487                  p=MINRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)                  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                  # solve Au=f-B^*p
1489                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1490                  #       u=v+(u-v)                  #       u=v+(u-v)
1491          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1492                else:
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 "DIFF=",util.integrate(p)
1518
1519                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1520
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