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

revision 1514 by artak, Wed Apr 16 00:15:44 2008 UTC revision 1552 by gross, Thu May 8 08:52:41 2008 UTC
# Line 477  class IterationHistory(object): Line 477  class IterationHistory(object):
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    ##############################################
585  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  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=20):  def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
598     iter=0     iter=0
# Line 597  def GMRESm(b, Aprod, Msolve, bilinearfor Line 602  def GMRESm(b, Aprod, Msolve, bilinearfor
602     norm_b=math.sqrt(r_dot_r)     norm_b=math.sqrt(r_dot_r)
603
604     if x==None:     if x==None:
605        x=0        x=0*b
606     else:     else:
607        r=Msolve(b-Aprod(x))        r=Msolve(b-Aprod(x))
608        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
# Line 613  def GMRESm(b, Aprod, Msolve, bilinearfor Line 618  def GMRESm(b, Aprod, Msolve, bilinearfor
618
619     v.append(r/rho)     v.append(r/rho)
620     g[0]=rho     g[0]=rho

621     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
622
623      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
# Line 628  def GMRESm(b, Aprod, Msolve, bilinearfor Line 632  def GMRESm(b, Aprod, Msolve, bilinearfor
632  # Modified Gram-Schmidt  # Modified Gram-Schmidt
633      for j in range(iter+1):      for j in range(iter+1):
634        h[j][iter]=bilinearform(v[j],v[iter+1])          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
641    # Reorthogonalize if needed
642        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
643         for j in range(iter+1):
644            hr=bilinearform(v[j],v[iter+1])
645                h[j][iter]=h[j][iter]+hr #vhat
646                v[iter+1] -= hr*v[j]
647
648         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
649         h[iter+1][iter]=v_norm2
650
651    #   watch out for happy breakdown
652            if not v_norm2 == 0:
653             v[iter+1]=v[iter+1]/h[iter+1][iter]
654
655    #   Form and store the information for the new Givens rotation
656        if iter > 0 :
657            hhat=numarray.zeros(iter+1,numarray.Float64)
658            for i in range(iter+1) : hhat[i]=h[i][iter]
659            hhat=givapp(c[0:iter],s[0:iter],hhat);
660                for i in range(iter+1) : h[i][iter]=hhat[i]
661
662        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
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            rho=abs(g[iter+1])
672        iter+=1
673
674    # At this point either iter > iter_max or rho < tol.
675    # It's time to compute x and leave.
676
677       if iter > 0 :
678         y=numarray.zeros(iter,numarray.Float64)
679         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
680         if iter > 1 :
681            i=iter-2
682            while i>=0 :
683              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
684              i=i-1
685         xhat=v[iter-1]*y[iter-1]
686         for i in range(iter-1):
687        xhat += v[i]*y[i]
688       else : xhat=v[0]
689
690       x += xhat
691       if iter<iter_restart-1:
692          stopped=True
693       else:
694          stopped=False
695
696       return x,stopped
697
698
699    ######################################################
700    def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
701    ######################################################
702
703    # DIRDER estimates the directional derivative of a function F.
704
705
706    # Hardwired difference increment.
707    #
708      epsnew = 1.0e-07
709    #
710    #  Scale the step.
711    #
712      norm_w=math.sqrt(bilinearform(w,w))
713      if norm_w== 0.0:
714        return x/x
715
716      epsnew = epsnew / norm_w
717
718      if norm_w > 0.0:
719        epsnew = epsnew * math.sqrt(bilinearform(x,x))
720    #
721    #  DEL and F1 could share the same space if storage
722    #  is more important than clarity.
723    #
724
725      DEL = x + epsnew * w
726      f1 = -Msolve(Aprod(DEL))
727      z = ( f1 - f0 ) / epsnew
728      return z
729
730    ######################################################
731    def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
732    ######################################################
733       b=-f0
734       b_dot_b = bilinearform(b, b)
735       if b_dot_b<0: raise NegativeNorm,"negative norm."
736       norm_b=math.sqrt(b_dot_b)
737
738       r=b
739
740       if x==None:
741          x=0*f0
742       else:
743          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0
744
745       r_dot_r = bilinearform(r, r)
746       if r_dot_r<0: raise NegativeNorm,"negative norm."
747
748       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
749       c=numarray.zeros(iter_restart,numarray.Float64)
750       s=numarray.zeros(iter_restart,numarray.Float64)
751       g=numarray.zeros(iter_restart,numarray.Float64)
752       v=[]
753
754       rho=math.sqrt(r_dot_r)
755
756       v.append(r/rho)
757       g[0]=rho
758       iter=0
759
760       while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
761
762        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
763
764
765            p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
766
767        v.append(p)
768
769        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
770
771    # Modified Gram-Schmidt
772        for j in range(iter+1):
773          h[j][iter]=bilinearform(v[j],v[iter+1])
774        v[iter+1]+=(-1.)*h[j][iter]*v[j]        v[iter+1]+=(-1.)*h[j][iter]*v[j]
775
776      h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))      h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
# Line 667  def GMRESm(b, Aprod, Msolve, bilinearfor Line 810  def GMRESm(b, Aprod, Msolve, bilinearfor
810          rho=abs(g[iter+1])          rho=abs(g[iter+1])
811      iter+=1      iter+=1
812
# At this point either iter > iter_max or rho < tol.
# It's time to compute x and leave.

813     if iter > 0 :     if iter > 0 :
814       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)
815       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1][iter-1]
# Line 691  def GMRESm(b, Aprod, Msolve, bilinearfor Line 831  def GMRESm(b, Aprod, Msolve, bilinearfor
831
832     return x,stopped     return x,stopped
833
834    #################################################
835  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
836    #################################################
837      #      #
838      #  minres solves the system of linear equations Ax = b      #  minres solves the system of linear equations Ax = b
839      #  where A is a symmetric matrix (possibly indefinite or singular)      #  where A is a symmetric matrix (possibly indefinite or singular)
# Line 722  def MINRES(b, Aprod, Msolve, bilinearfor Line 863  def MINRES(b, Aprod, Msolve, bilinearfor
863
864      r1    = b      r1    = b
865      y = Msolve(b)      y = Msolve(b)
866      beta1 = bilinearform(b,y)      beta1 = bilinearform(y,b)
867
868      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
869
# Line 758  def MINRES(b, Aprod, Msolve, bilinearfor Line 899  def MINRES(b, Aprod, Msolve, bilinearfor
899      #---------------------------------------------------------------------      #---------------------------------------------------------------------
900      # Main iteration loop.      # Main iteration loop.
901      # --------------------------------------------------------------------      # --------------------------------------------------------------------
902      while not stoppingcriterium(rnorm,Anorm*ynorm):    #  checks ||r|| < (||A|| ||x||) * TOL      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
903
904      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
905          iter    = iter  +  1          iter    = iter  +  1
# Line 784  def MINRES(b, Aprod, Msolve, bilinearfor Line 925  def MINRES(b, Aprod, Msolve, bilinearfor
925            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
926
927          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
928          y      = (- alfa/beta)*r2 + y          y      += (- alfa/beta)*r2
929          r1     = r2          r1     = r2
930          r2     = y          r2     = y
931          y = Msolve(r2)          y = Msolve(r2)
932          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
933          beta   = bilinearform(r2,y)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
934          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm,"negative norm."
935
936          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
# Line 824  def MINRES(b, Aprod, Msolve, bilinearfor Line 965  def MINRES(b, Aprod, Msolve, bilinearfor
965          w1    = w2          w1    = w2
966          w2    = w          w2    = w
967          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
968          x     = x  +  phi*w          x     +=  phi*w
969
970          # Go round again.          # Go round again.
971
# Line 843  def MINRES(b, Aprod, Msolve, bilinearfor Line 984  def MINRES(b, Aprod, Msolve, bilinearfor
984          rnorm  = phibar          rnorm  = phibar
985
986      return x      return x
987
988    ######################################
989    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
990    #####################################
991        gamma=.9
992        lmaxit=100
993        etamax=.5
994
995        n = 1 #len(x)
996        iter=0
997
998        # evaluate f at the initial iterate
999        # compute the stop tolerance
1000        #
1001        r=b
1002        if x==None:
1003          x=0*b
1004        else:
1005          r += (-1)*Aprod(x)
1006
1007        f0=-Msolve(r)
1008        fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1009        fnrmo=1
1010        stop_tol=atol + rtol*fnrm
1011        #
1012        # main iteration loop
1013        #
1014        while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):
1015
1016                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1017            #
1018            # keep track of the ratio (rat = fnrm/frnmo)
1019            # of successive residual norms and
1020            # the iteration counter (iter)
1021            #
1022            #rat=fnrm/fnrmo
1023            fnrmo=fnrm
1024            iter+=1
1025            #
1026                # compute the step using a GMRES(m) routine especially designed for this purpose
1027            #
1028                initer=0
1029                while True:
1030                   xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1031                   if stopped: break
1032                   initer+=iter_restart
1033            x+=xc
1034            f0=-Msolve(Aprod(x))
1035            fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1036            rat=fnrm/fnrmo
1037
1038
1040        #
1041            if etamax > 0:
1042                etaold=etamax
1043                etanew=gamma*rat*rat
1044                if gamma*etaold*etaold > .1 :
1045                    etanew=max(etanew,gamma*etaold*etaold)
1046
1047                etamax=min(etanew,etamax)
1048                etamax=max(etamax,.5*stop_tol/fnrm)
1049
1050        return x
1051
1052  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1053
# Line 888  def TFQMR(b, Aprod, Msolve, bilinearform Line 1092  def TFQMR(b, Aprod, Msolve, bilinearform
1092  #  #
1093  #  while ( iter < kmax-1 ):  #  while ( iter < kmax-1 ):
1094
1095    while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b):    while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1096      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
1097
1098      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
# Line 944  def TFQMR(b, Aprod, Msolve, bilinearform Line 1148  def TFQMR(b, Aprod, Msolve, bilinearform
1148
1149      iter = iter + 1      iter = iter + 1
1150
print iter
1151    return x    return x
1152
1153
# Line 1033  class ArithmeticTuple(object): Line 1236  class ArithmeticTuple(object):
1236         """         """
1237         out=[]         out=[]
1238         for i in range(len(self)):         for i in range(len(self)):
1239             out.append(self[i]/other)             out.append(self[i]*(1./other))
1240         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1241
1242     def __rdiv__(self,other):     def __rdiv__(self,other):
# Line 1065  class ArithmeticTuple(object): Line 1268  class ArithmeticTuple(object):
1268             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1269         return self         return self
1270
1272           """
1274
1275           @param other: increment
1276           @type other: C{ArithmeticTuple}
1277           """
1278           if not isinstance(other,float):
1279            if len(self) != len(other):
1280               raise ValueError,"tuple length must match."
1281            for i in range(len(self)):
1282               self.__items[i]+=other[i]
1283           else:
1284            for i in range(len(self)):
1285               self.__items[i]+=other
1286
1287           return self
1288
1289       def __sub__(self,other):
1290           """
1291           subtract other from self
1292
1293           @param other: increment
1294           @type other: C{ArithmeticTuple}
1295           """
1296           if len(self) != len(other):
1297               raise ValueError,"tuple length must match."
1298           for i in range(len(self)):
1299               self.__items[i]-=other[i]
1300           return self
1301
1302       def __neg__(self):
1303           """
1304           negate
1305
1306           """
1307           for i in range(len(self)):
1308               self.__items[i]=-self.__items[i]
1309           return self
1310
1311
1313        """        """
1314        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
1380
1381        def __inner_p(self,p1,p2):        def __inner_p(self,p1,p2):
1382           return self.inner(p1,p2)           return self.inner(p1,p2)
1383
1384          def __inner_a(self,a1,a2):
1385             return self.inner_a(a1,a2)
1386
1387          def __inner_a1(self,a1,a2):
1388             return self.inner(a1[1],a2[1])
1389
1390        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
1391            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
1392
1393        def __stoppingcriterium_GMRES(self,norm_r,norm_b):        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1394            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)

1395
1396        def setTolerance(self,tolerance=1.e-8):        def setTolerance(self,tolerance=1.e-8):
1397                self.__tol=tolerance                self.__tol=tolerance
1402        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1403                return self.__reduction*self.getTolerance()                return self.__reduction*self.getTolerance()
1404
1405        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=10):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1406                """                """
1407                solves the saddle point problem using initial guesses v and p.                solves the saddle point problem using initial guesses v and p.
1408
1415                # which leads to BA^-1B^*p = BA^-1f                  # which leads to BA^-1B^*p = BA^-1f
1416
1417            # 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)

1418            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)

1419                Bz=self.B(self.__z)                Bz=self.B(self.__z)
1420                #                #
1421            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1422                #                #
#   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1423                #                #
#   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)
1424                #                #
1425                self.iter=0                self.iter=0
1426            if solver=='GMRES':                  if solver=='GMRES':
1427                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1428                  p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)                  p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1429                    # solve Au=f-B^*p
1430                    #       A(u-v)=f-B^*p-Av
1431                    #       u=v+(u-v)
1432            u=v+self.solve_A(v,p)
1433
1434              if solver=='NewtonGMRES':
1435                    if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1436                    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())
1437                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1438                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1439                  #       u=v+(u-v)                  #       u=v+(u-v)
1440          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1441
1442
1443            if solver=='TFQMR':                  if solver=='TFQMR':
1444                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1445                  p=TFQMR(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)                  p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1446                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1447                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1448                  #       u=v+(u-v)                  #       u=v+(u-v)
1450
1451            if solver=='MINRES':                  if solver=='MINRES':
1452                  if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter                  if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1453                  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.)
1454                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1455                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1456                  #       u=v+(u-v)                  #       u=v+(u-v)
1457          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1458
1459              if solver=='GMRESC':
1460                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1461                    p0=self.solve_prec(Bz)
1462               #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1463                   #p-=alfa
1464                    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)
1465                    #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())
1466                    # solve Au=f-B^*p
1467                    #       A(u-v)=f-B^*p-Av
1468                    #       u=v+(u-v)
1469                p=x[1]
1470            #u=v+self.solve_A(v,p)
1471            #p=x[1]
1472            #u=x[0]
1473
1474                if solver=='PCG':                if solver=='PCG':
1475                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1476                    #
1477                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1478                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1479                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1480                  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)
1481              u=r[0]                u=r[0]
1482                    print "DIFF=",util.integrate(p)
1483
1484                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1485
1488        def __Msolve(self,r):        def __Msolve(self,r):
1489            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1490
1491        def __Msolve_GMRES(self,r):        def __Msolve2(self,r):
1492            return self.solve_prec(r)            return self.solve_prec(r*1)
1493
1494          def __Mempty(self,r):
1495              return r
1496
1497
1498        def __Aprod(self,p):        def __Aprod(self,p):
1499            # return BA^-1B*p            # return BA^-1B*p
1500            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1501            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1502            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1503
1504        def __Aprod_GMRES(self,p):        def __Anext(self,x):
1505            # return BA^-1B*p            # return next v,p
1506            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
1507
1508          pc=x[1]
1509              v=self.solve_A(self.__z,-pc)
1510          p=self.solve_prec(self.B(v))
1511
1512              return ArithmeticTuple(v,p)
1513
1514
1515          def __Aprod2(self,p):
1516              # return BA^-1B*p
1517              #solve Av =B^*p as Av =f-Az-B^*(-p)
1518        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1519            return self.B(v)            return self.B(v)
1520
1521          def __Aprod_Newton(self,p):
1522              # return BA^-1B*p - Bz
1523              #solve Av =-B^*p as Av =f-Az-B^*p
1524          v=self.solve_A(self.__z,-p)
1525              return self.B(v-self.__z)
1526
1527          def __Aprod_Newton2(self,x):
1528              # return BA^-1B*p - Bz
1529              #solve Av =-B^*p as Av =f-Az-B^*p
1530              pc=x[1]
1531          v=self.solve_A(self.__z,-pc)
1532              p=self.solve_prec(self.B(v-self.__z))
1533              return ArithmeticTuple(v,p)
1534