/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1414 by gross, Thu Feb 14 10:01:43 2008 UTC revision 1519 by artak, Tue Apr 22 03:45:36 2008 UTC
# Line 48  import numarray Line 48  import numarray
48  import util  import util
49  import math  import math
50    
51    ##### Added by Artak
52    # from Numeric import zeros,Int,Float64
53    ###################################
54    
55    
56  class TimeIntegrationManager:  class TimeIntegrationManager:
57    """    """
58    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
# Line 469  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,solver="GMRES",TOL=None):
481           """
482           returns True if the C{norm_r} is C{tolerance}*C{norm_b}
483    
484          
485           @param norm_r: current residual norm
486           @type norm_r: non-negative C{float}
487           @param norm_b: norm of right hand side
488           @type norm_b: non-negative C{float}
489           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
490           @rtype: C{bool}
491    
492           """
493           if TOL==None:
494              TOL=self.tolerance
495           self.history.append(norm_r)
496           if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])
497           return self.history[-1]<=TOL * norm_b
498    
499  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
500     """     """
501     Solver for     Solver for
# Line 539  type like argument C{x}. Line 563  type like argument C{x}.
563    
564     return x,r     return x,r
565    
566    
567    ############################
568    # Added by Artak
569    #################################3
570    
571    #Apply a sequence of k Givens rotations, used within gmres codes
572    # vrot=givapp(c, s, vin, k)
573    def givapp(c,s,vin):
574        vrot=vin # warning: vin is altered!!!!
575        if isinstance(c,float):
576            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
577        else:
578            for i in range(len(c)):
579                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
580            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
581                vrot[i:i+2]=w1,w2
582        return vrot
583    
584    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
585       m=iter_restart
586       iter=0
587       while True:
588          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
589          x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
590          if stopped: break
591          iter+=iter_restart    
592       return x
593    
594    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
595       iter=0
596       r=Msolve(b)
597       r_dot_r = bilinearform(r, r)
598       if r_dot_r<0: raise NegativeNorm,"negative norm."
599       norm_b=math.sqrt(r_dot_r)
600    
601       if x==None:
602          x=0
603       else:
604          r=Msolve(b-Aprod(x))
605          r_dot_r = bilinearform(r, r)
606          if r_dot_r<0: raise NegativeNorm,"negative norm."
607      
608       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
609       c=numarray.zeros(iter_restart,numarray.Float64)
610       s=numarray.zeros(iter_restart,numarray.Float64)
611       g=numarray.zeros(iter_restart,numarray.Float64)
612       v=[]
613    
614       rho=math.sqrt(r_dot_r)
615      
616       v.append(r/rho)
617       g[0]=rho
618       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
619    
620        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
621    
622        
623        p=Msolve(Aprod(v[iter]))
624    
625        v.append(p)
626    
627        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
628    
629    # Modified Gram-Schmidt
630        for j in range(iter+1):
631          h[j][iter]=bilinearform(v[j],v[iter+1])  
632          v[iter+1]+=(-1.)*h[j][iter]*v[j]
633          
634        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
635        v_norm2=h[iter+1][iter]
636    
637    
638    # Reorthogonalize if needed
639        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
640         for j in range(iter+1):
641            hr=bilinearform(v[j],v[iter+1])
642                h[j][iter]=h[j][iter]+hr #vhat
643                v[iter+1] +=(-1.)*hr*v[j]
644    
645         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
646         h[iter+1][iter]=v_norm2
647    
648    #   watch out for happy breakdown
649            if v_norm2 != 0:
650             v[iter+1]=v[iter+1]/h[iter+1][iter]
651    
652    #   Form and store the information for the new Givens rotation
653        if iter > 0 :
654            hhat=[]
655            for i in range(iter+1) : hhat.append(h[i][iter])
656            hhat=givapp(c[0:iter],s[0:iter],hhat);
657                for i in range(iter+1) : h[i][iter]=hhat[i]
658    
659        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
660        if mu!=0 :
661            c[iter]=h[iter][iter]/mu
662            s[iter]=-h[iter+1][iter]/mu
663            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
664            h[iter+1][iter]=0.0
665            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
666    
667    # Update the residual norm
668            rho=abs(g[iter+1])
669        iter+=1
670    
671    # At this point either iter > iter_max or rho < tol.
672    # It's time to compute x and leave.        
673    
674       if iter > 0 :
675         y=numarray.zeros(iter,numarray.Float64)    
676         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
677         if iter > 1 :  
678            i=iter-2  
679            while i>=0 :
680              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
681              i=i-1
682         xhat=v[iter-1]*y[iter-1]
683         for i in range(iter-1):
684        xhat += v[i]*y[i]
685       else : xhat=v[0]
686        
687       x += xhat
688       if iter<iter_restart-1:
689          stopped=True
690       else:
691          stopped=False
692    
693       return x,stopped
694    
695    ######################################################
696    def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
697    ######################################################
698    
699    # DIRDER estimates the directional derivative of a function F.
700    
701    
702    # Hardwired difference increment.
703    #
704      epsnew = 1.0e-07
705    #
706    #  Scale the step.
707    #
708      norm_w=math.sqrt(bilinearform(w,w))
709      if norm_w== 0.0:
710        return x/x
711    
712      epsnew = epsnew / norm_w
713    
714      if norm_w > 0.0:
715        epsnew = epsnew * math.sqrt(bilinearform(x,x))
716    #
717    #  DEL and F1 could share the same space if storage
718    #  is more important than clarity.
719    #
720    
721      DEL = x + epsnew * w
722      f1 = -Msolve(Aprod(DEL))
723      z = ( f1 - f0 ) / epsnew
724      return z
725    
726    ######################################################
727    def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
728    ######################################################
729       b=-f0
730       b_dot_b = bilinearform(b, b)
731       if b_dot_b<0: raise NegativeNorm,"negative norm."
732       norm_b=math.sqrt(b_dot_b)
733    
734       r=b
735    
736       if x==None:
737          x=0
738       else:
739          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0  
740          
741       r_dot_r = bilinearform(r, r)
742       if r_dot_r<0: raise NegativeNorm,"negative norm."
743      
744       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
745       c=numarray.zeros(iter_restart,numarray.Float64)
746       s=numarray.zeros(iter_restart,numarray.Float64)
747       g=numarray.zeros(iter_restart,numarray.Float64)
748       v=[]
749    
750       rho=math.sqrt(r_dot_r)
751      
752       v.append(r/rho)
753       g[0]=rho
754       iter=0
755    
756       while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
757    
758        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
759    
760        
761            p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
762    
763        v.append(p)
764    
765        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
766    
767    # Modified Gram-Schmidt
768        for j in range(iter+1):
769          h[j][iter]=bilinearform(v[j],v[iter+1])  
770          v[iter+1]+=(-1.)*h[j][iter]*v[j]
771          
772        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
773        v_norm2=h[iter+1][iter]
774    
775    
776    # Reorthogonalize if needed
777        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
778         for j in range(iter+1):
779            hr=bilinearform(v[j],v[iter+1])
780                h[j][iter]=h[j][iter]+hr #vhat
781                v[iter+1] +=(-1.)*hr*v[j]
782    
783         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
784         h[iter+1][iter]=v_norm2
785    
786    #   watch out for happy breakdown
787            if v_norm2 != 0:
788             v[iter+1]=v[iter+1]/h[iter+1][iter]
789    
790    #   Form and store the information for the new Givens rotation
791        if iter > 0 :
792            hhat=[]
793            for i in range(iter+1) : hhat.append(h[i][iter])
794            hhat=givapp(c[0:iter],s[0:iter],hhat);
795                for i in range(iter+1) : h[i][iter]=hhat[i]
796    
797        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
798        if mu!=0 :
799            c[iter]=h[iter][iter]/mu
800            s[iter]=-h[iter+1][iter]/mu
801            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
802            h[iter+1][iter]=0.0
803            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
804    
805    # Update the residual norm
806            rho=abs(g[iter+1])
807        iter+=1
808    
809       if iter > 0 :
810         y=numarray.zeros(iter,numarray.Float64)    
811         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
812         if iter > 1 :  
813            i=iter-2  
814            while i>=0 :
815              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
816              i=i-1
817         xhat=v[iter-1]*y[iter-1]
818         for i in range(iter-1):
819        xhat += v[i]*y[i]
820       else : xhat=v[0]
821        
822       x += xhat
823       if iter<iter_restart-1:
824          stopped=True
825       else:
826          stopped=False
827    
828       return x,stopped
829    
830    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
831    
832        #
833        #  minres solves the system of linear equations Ax = b
834        #  where A is a symmetric matrix (possibly indefinite or singular)
835        #  and b is a given vector.
836        #  
837        #  "A" may be a dense or sparse matrix (preferably sparse!)
838        #  or the name of a function such that
839        #               y = A(x)
840        #  returns the product y = Ax for any given vector x.
841        #
842        #  "M" defines a positive-definite preconditioner M = C C'.
843        #  "M" may be a dense or sparse matrix (preferably sparse!)
844        #  or the name of a function such that
845        #  solves the system My = x for any given vector x.
846        #
847        #
848        
849        #------------------------------------------------------------------
850        # Set up y and v for the first Lanczos vector v1.
851        # y  =  beta1 P' v1,  where  P = C**(-1).
852        # v is really P' v1.
853        #------------------------------------------------------------------
854        if x==None:
855          x=0*b
856        else:
857          b += (-1)*Aprod(x)
858    
859        r1    = b
860        y = Msolve(b)
861        beta1 = bilinearform(b,y)
862    
863        if beta1< 0: raise NegativeNorm,"negative norm."
864    
865        #  If b = 0 exactly, stop with x = 0.
866        if beta1==0: return x*0.
867    
868        if beta1> 0:
869          beta1  = math.sqrt(beta1)      
870    
871        #------------------------------------------------------------------
872        # Initialize quantities.
873        # ------------------------------------------------------------------
874        iter   = 0
875        Anorm = 0
876        ynorm = 0
877        oldb   = 0
878        beta   = beta1
879        dbar   = 0
880        epsln  = 0
881        phibar = beta1
882        rhs1   = beta1
883        rhs2   = 0
884        rnorm  = phibar
885        tnorm2 = 0
886        ynorm2 = 0
887        cs     = -1
888        sn     = 0
889        w      = b*0.
890        w2     = b*0.
891        r2     = r1
892        eps    = 0.0001
893    
894        #---------------------------------------------------------------------
895        # Main iteration loop.
896        # --------------------------------------------------------------------
897        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
898    
899        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
900            iter    = iter  +  1
901    
902            #-----------------------------------------------------------------
903            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
904            # The general iteration is similar to the case k = 1 with v0 = 0:
905            #
906            #   p1      = Operator * v1  -  beta1 * v0,
907            #   alpha1  = v1'p1,
908            #   q2      = p2  -  alpha1 * v1,
909            #   beta2^2 = q2'q2,
910            #   v2      = (1/beta2) q2.
911            #
912            # Again, y = betak P vk,  where  P = C**(-1).
913            #-----------------------------------------------------------------
914            s = 1/beta                 # Normalize previous vector (in y).
915            v = s*y                    # v = vk if P = I
916        
917            y      = Aprod(v)
918        
919            if iter >= 2:
920              y = y - (beta/oldb)*r1
921    
922            alfa   = bilinearform(v,y)              # alphak
923            y      = (- alfa/beta)*r2 + y
924            r1     = r2
925            r2     = y
926            y = Msolve(r2)
927            oldb   = beta                           # oldb = betak
928            beta   = bilinearform(r2,y)             # beta = betak+1^2
929            if beta < 0: raise NegativeNorm,"negative norm."
930    
931            beta   = math.sqrt( beta )
932            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
933            
934            if iter==1:                 # Initialize a few things.
935              gmax   = abs( alfa )      # alpha1
936              gmin   = gmax             # alpha1
937    
938            # Apply previous rotation Qk-1 to get
939            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
940            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
941        
942            oldeps = epsln
943            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
944            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
945            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
946            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
947    
948            # Compute the next plane rotation Qk
949    
950            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
951            gamma  = max(gamma,eps)
952            cs     = gbar / gamma             # ck
953            sn     = beta / gamma             # sk
954            phi    = cs * phibar              # phik
955            phibar = sn * phibar              # phibark+1
956    
957            # Update  x.
958    
959            denom = 1/gamma
960            w1    = w2
961            w2    = w
962            w     = (v - oldeps*w1 - delta*w2) * denom
963            x     = x  +  phi*w
964    
965            # Go round again.
966    
967            gmax   = max(gmax,gamma)
968            gmin   = min(gmin,gamma)
969            z      = rhs1 / gamma
970            ynorm2 = z*z  +  ynorm2
971            rhs1   = rhs2 -  delta*z
972            rhs2   =      -  epsln*z
973    
974            # Estimate various norms and test for convergence.
975    
976            Anorm  = math.sqrt( tnorm2 )
977            ynorm  = math.sqrt( ynorm2 )
978    
979            rnorm  = phibar
980    
981        return x
982        
983    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20):
984    
985        gamma=.9
986        lmaxit=40
987        etamax=.5
988    
989        n = 1 #len(x)
990        iter=0
991        
992        # evaluate f at the initial iterate
993        # compute the stop tolerance
994        #
995        r=b
996        if x==None:
997          x=0*b
998        else:
999          r += (-1)*Aprod(x)
1000    
1001        f0=-Msolve(r)
1002        fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1003        fnrmo=1
1004        atol=1.e-2
1005        rtol=1.e-4
1006        stop_tol=atol + rtol*fnrm
1007        #
1008        # main iteration loop
1009        #
1010        while not stoppingcriterium(fnrm,stop_tol,'NewtonGMRES',TOL=1.):
1011    
1012                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1013            #
1014            # keep track of the ratio (rat = fnrm/frnmo)
1015            # of successive residual norms and
1016            # the iteration counter (iter)
1017            #
1018            #rat=fnrm/fnrmo
1019            fnrmo=fnrm
1020            iter+=1
1021            #
1022                # compute the step using a GMRES(m) routine especially designed for this purpose
1023            #
1024                initer=0
1025                while True:
1026                   xc,stopped=FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1027                   if stopped: break
1028                   initer+=iter_restart
1029            xold=x
1030            x+=xc
1031            f0=-Msolve(Aprod(x))
1032            fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1033            rat=fnrm/fnrmo
1034    
1035    
1036        #   adjust eta
1037        #
1038            if etamax > 0:
1039                etaold=etamax
1040                etanew=gamma*rat*rat
1041                if gamma*etaold*etaold > .1 :
1042                    etanew=max(etanew,gamma*etaold*etaold)
1043                
1044                etamax=min(etanew,etamax)
1045                etamax=max(etamax,.5*stop_tol/fnrm)
1046    
1047        return x
1048    
1049    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1050    
1051    # TFQMR solver for linear systems
1052    #
1053    #
1054    # initialization
1055    #
1056      errtol = math.sqrt(bilinearform(b,b))
1057      norm_b=errtol
1058      kmax  = iter_max
1059      error = []
1060    
1061      if math.sqrt(bilinearform(x,x)) != 0.0:
1062        r = b - Aprod(x)
1063      else:
1064        r = b
1065    
1066      r=Msolve(r)
1067    
1068      u1=0
1069      u2=0
1070      y1=0
1071      y2=0
1072    
1073      w = r
1074      y1 = r
1075      iter = 0
1076      d = 0
1077      
1078      v = Msolve(Aprod(y1))
1079      u1 = v
1080      
1081      theta = 0.0;
1082      eta = 0.0;
1083      tau = math.sqrt(bilinearform(r,r))
1084      error = [ error, tau ]
1085      rho = tau * tau
1086      m=1
1087    #
1088    #  TFQMR iteration
1089    #
1090    #  while ( iter < kmax-1 ):
1091      
1092      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1093        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1094    
1095        sigma = bilinearform(r,v)
1096    
1097        if ( sigma == 0.0 ):
1098          raise 'TFQMR breakdown, sigma=0'
1099        
1100    
1101        alpha = rho / sigma
1102    
1103        for j in range(2):
1104    #
1105    #   Compute y2 and u2 only if you have to
1106    #
1107          if ( j == 1 ):
1108            y2 = y1 - alpha * v
1109            u2 = Msolve(Aprod(y2))
1110    
1111          m = 2 * (iter+1) - 2 + (j+1)
1112          if j==0:
1113             w = w - alpha * u1
1114             d = y1 + ( theta * theta * eta / alpha ) * d
1115          if j==1:
1116             w = w - alpha * u2
1117             d = y2 + ( theta * theta * eta / alpha ) * d
1118    
1119          theta = math.sqrt(bilinearform(w,w))/ tau
1120          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1121          tau = tau * theta * c
1122          eta = c * c * alpha
1123          x = x + eta * d
1124    #
1125    #  Try to terminate the iteration at each pass through the loop
1126    #
1127         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1128         #   error = [ error, tau ]
1129         #   total_iters = iter
1130         #   break
1131          
1132    
1133        if ( rho == 0.0 ):
1134          raise 'TFQMR breakdown, rho=0'
1135        
1136    
1137        rhon = bilinearform(r,w)
1138        beta = rhon / rho;
1139        rho = rhon;
1140        y1 = w + beta * y2;
1141        u1 = Msolve(Aprod(y1))
1142        v = u1 + beta * ( u2 + beta * v )
1143        error = [ error, tau ]
1144        total_iters = iter
1145        
1146        iter = iter + 1
1147    
1148      return x
1149    
1150    
1151    #############################################
1152    
1153  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1154     """     """
1155     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
# Line 608  class ArithmeticTuple(object): Line 1219  class ArithmeticTuple(object):
1219             out.append(other*self[i])             out.append(other*self[i])
1220         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1221    
1222    #########################
1223    # Added by Artak
1224    #########################
1225       def __div__(self,other):
1226           """
1227           dividing from the right
1228    
1229           @param other: scaling factor
1230           @type other: C{float}
1231           @return: itemwise self/other
1232           @rtype: L{ArithmeticTuple}
1233           """
1234           out=[]
1235           for i in range(len(self)):
1236               out.append(self[i]/other)
1237           return ArithmeticTuple(*tuple(out))
1238    
1239       def __rdiv__(self,other):
1240           """
1241           dividing from the left
1242    
1243           @param other: scaling factor
1244           @type other: C{float}
1245           @return: itemwise other/self
1246           @rtype: L{ArithmeticTuple}
1247           """
1248           out=[]
1249           for i in range(len(self)):
1250               out.append(other/self[i])
1251           return ArithmeticTuple(*tuple(out))
1252      
1253    ##########################################33
1254    
1255     def __iadd__(self,other):     def __iadd__(self,other):
1256         """         """
1257         in-place add of other to self         in-place add of other to self
# Line 690  class HomogeneousSaddlePointProblem(obje Line 1334  class HomogeneousSaddlePointProblem(obje
1334        def __inner(self,p,r):        def __inner(self,p,r):
1335           return self.inner(p,r[1])           return self.inner(p,r[1])
1336    
1337          def __inner_p(self,p1,p2):
1338             return self.inner(p1,p2)
1339    
1340        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
1341            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
1342    
1343          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1344              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1345    
1346        def setTolerance(self,tolerance=1.e-8):        def setTolerance(self,tolerance=1.e-8):
1347                self.__tol=tolerance                self.__tol=tolerance
1348        def getTolerance(self):        def getTolerance(self):
# Line 702  class HomogeneousSaddlePointProblem(obje Line 1352  class HomogeneousSaddlePointProblem(obje
1352        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1353                return self.__reduction*self.getTolerance()                return self.__reduction*self.getTolerance()
1354    
1355        def solve(self,v,p,max_iter=20, verbose=False, show_details=False):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1356                """                """
1357                solves the saddle point problem using initial guesses v and p.                solves the saddle point problem using initial guesses v and p.
1358    
# Line 711  class HomogeneousSaddlePointProblem(obje Line 1361  class HomogeneousSaddlePointProblem(obje
1361                self.verbose=verbose                self.verbose=verbose
1362                self.show_details=show_details and self.verbose                self.show_details=show_details and self.verbose
1363    
1364                  # assume p is known: then v=A^-1(f-B^*p)
1365                  # which leads to BA^-1B^*p = BA^-1f  
1366    
1367            # 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)      
1368    
1369              
1370            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
1371    
1372                Bz=self.B(self.__z)                Bz=self.B(self.__z)
1373                #                #
1374            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
# Line 724  class HomogeneousSaddlePointProblem(obje Line 1379  class HomogeneousSaddlePointProblem(obje
1379                #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)                #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1380                #                #
1381                self.iter=0                self.iter=0
1382                if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter            if solver=='GMRES':      
1383                p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1384            return r[0],p                  p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1385                    # solve Au=f-B^*p
1386                    #       A(u-v)=f-B^*p-Av
1387                    #       u=v+(u-v)
1388            u=v+self.solve_A(v,p)
1389    
1390              if solver=='NewtonGMRES':    
1391                    if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1392                    p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1393                    # solve Au=f-B^*p
1394                    #       A(u-v)=f-B^*p-Av
1395                    #       u=v+(u-v)
1396            u=v+self.solve_A(v,p)
1397                    
1398    
1399              if solver=='TFQMR':      
1400                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1401                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1402                    # solve Au=f-B^*p
1403                    #       A(u-v)=f-B^*p-Av
1404                    #       u=v+(u-v)
1405            u=v+self.solve_A(v,p)
1406    
1407              if solver=='MINRES':      
1408                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1409                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1410                    # solve Au=f-B^*p
1411                    #       A(u-v)=f-B^*p-Av
1412                    #       u=v+(u-v)
1413            u=v+self.solve_A(v,p)
1414                  
1415                  if solver=='PCG':
1416                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1417                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1418                u=r[0]  
1419    
1420                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1421    
1422              return u,p
1423    
1424        def __Msolve(self,r):        def __Msolve(self,r):
1425            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1426    
1427          def __Msolve2(self,r):
1428              return self.solve_prec(r)
1429    
1430    
1431        def __Aprod(self,p):        def __Aprod(self,p):
1432            # return BA^-1B*p            # return BA^-1B*p
1433            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
1434            v=self.solve_A(self.__z,p)            v=self.solve_A(self.__z,-p)
1435            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1436    
1437          def __Aprod2(self,p):
1438              # return BA^-1B*p
1439              #solve Av =-B^*p as Av =f-Az-B^*p
1440          v=self.solve_A(self.__z,-p)
1441              return self.B(v)
1442    
1443          def __Aprod_Newton(self,p):
1444              # return BA^-1B*p
1445              #solve Av =-B^*p as Av =f-Az-B^*p
1446          v=self.solve_A(self.__z,-p)
1447              return self.B(v-self.__z)
1448    
1449  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1450     """     """
# Line 775  class SaddlePointProblem(object): Line 1483  class SaddlePointProblem(object):
1483         @param text: a text message         @param text: a text message
1484         @type text: C{str}         @type text: C{str}
1485         """         """
1486         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "#s: #s"%(str(self),text)
1487    
1488     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1489         """         """
# Line 993  def MaskFromBoundaryTag(function_space,* Line 1701  def MaskFromBoundaryTag(function_space,*
1701        return util.whereNonZero(util.interpolate(out,function_space))        return util.whereNonZero(util.interpolate(out,function_space))
1702    
1703    
1704    

Legend:
Removed from v.1414  
changed lines
  Added in v.1519

  ViewVC Help
Powered by ViewVC 1.1.26