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

revision 1480 by gross, Mon Apr 7 23:38:50 2008 UTC revision 1481 by artak, Wed Apr 9 00:45:47 2008 UTC
# Line 689  def GMRESm(b, Aprod, Msolve, bilinearfor Line 689  def GMRESm(b, Aprod, Msolve, bilinearfor
689        stopped=False        stopped=False
690
691     return x,stopped     return x,stopped
692
693    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
694
695        #
696        #  minres solves the system of linear equations Ax = b
697        #  where A is a symmetric matrix (possibly indefinite or singular)
698        #  and b is a given vector.
699        #
700        #  "A" may be a dense or sparse matrix (preferably sparse!)
701        #  or the name of a function such that
702        #               y = A(x)
703        #  returns the product y = Ax for any given vector x.
704        #
705        #  "M" defines a positive-definite preconditioner M = C C'.
706        #  "M" may be a dense or sparse matrix (preferably sparse!)
707        #  or the name of a function such that
708        #  solves the system My = x for any given vector x.
709        #
710        #
711
712        #  Initialize
713
714        iter   = 0
715        Anorm = 0
716        ynorm = 0
717        x=x*0
718        #------------------------------------------------------------------
719        # Set up y and v for the first Lanczos vector v1.
720        # y  =  beta1 P' v1,  where  P = C**(-1).
721        # v is really P' v1.
722        #------------------------------------------------------------------
723        r1    = b
724        y = Msolve(b)
725        beta1 = bilinearform(b,y)
726
727        if beta1< 0: raise NegativeNorm,"negative norm."
728
729        #  If b = 0 exactly, stop with x = 0.
730        if beta1==0: return x*0.
731
732        if beta1> 0:
733          beta1  = math.sqrt(beta1)       # Normalize y to get v1 later.
734
735        #------------------------------------------------------------------
736        # Initialize other quantities.
737        # ------------------------------------------------------------------
738        oldb   = 0
739        beta   = beta1
740        dbar   = 0
741        epsln  = 0
742        phibar = beta1
743        rhs1   = beta1
744        rhs2   = 0
745        rnorm  = phibar
746        tnorm2 = 0
747        ynorm2 = 0
748        cs     = -1
749        sn     = 0
750        w      = b*0.
751        w2     = b*0.
752        r2     = r1
753        eps    = 0.0001
754
755        #---------------------------------------------------------------------
756        # Main iteration loop.
757        # --------------------------------------------------------------------
758        while not stoppingcriterium(rnorm,Anorm*ynorm):    #  ||r|| / (||A|| ||x||)
759
760        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
761            iter    = iter  +  1
762
763            #-----------------------------------------------------------------
764            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
765            # The general iteration is similar to the case k = 1 with v0 = 0:
766            #
767            #   p1      = Operator * v1  -  beta1 * v0,
768            #   alpha1  = v1'p1,
769            #   q2      = p2  -  alpha1 * v1,
770            #   beta2^2 = q2'q2,
771            #   v2      = (1/beta2) q2.
772            #
773            # Again, y = betak P vk,  where  P = C**(-1).
774            #-----------------------------------------------------------------
775            s = 1/beta                 # Normalize previous vector (in y).
776            v = s*y                    # v = vk if P = I
777
778            y      = Aprod(v)
779
780            if iter >= 2:
781              y = y - (beta/oldb)*r1
782
783            alfa   = bilinearform(v,y)              # alphak
784            y      = (- alfa/beta)*r2 + y
785            r1     = r2
786            r2     = y
787            y = Msolve(r2)
788            oldb   = beta                           # oldb = betak
789            beta   = bilinearform(r2,y)             # beta = betak+1^2
790            if beta < 0: raise NegativeNorm,"negative norm."
791
792            beta   = math.sqrt( beta )
793            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
794
795            if iter==1:                 # Initialize a few things.
796              gmax   = abs( alfa )      # alpha1
797              gmin   = gmax             # alpha1
798
799            # Apply previous rotation Qk-1 to get
800            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
801            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
802
803            oldeps = epsln
804            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
805            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
806            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
807            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
808
809            # Compute the next plane rotation Qk
810
811            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
812            gamma  = max(gamma,eps)
813            cs     = gbar / gamma             # ck
814            sn     = beta / gamma             # sk
815            phi    = cs * phibar              # phik
816            phibar = sn * phibar              # phibark+1
817
818            # Update  x.
819
820            denom = 1/gamma
821            w1    = w2
822            w2    = w
823            w     = (v - oldeps*w1 - delta*w2) * denom
824            x     = x  +  phi*w
825
826            # Go round again.
827
828            gmax   = max(gmax,gamma)
829            gmin   = min(gmin,gamma)
830            z      = rhs1 / gamma
831            ynorm2 = z*z  +  ynorm2
832            rhs1   = rhs2 -  delta*z
833            rhs2   =      -  epsln*z
834
835            # Estimate various norms and test for convergence.
836
837            Anorm  = math.sqrt( tnorm2 )
838            ynorm  = math.sqrt( ynorm2 )
839
840            rnorm  = phibar
841
843        return x
844
845  #############################################  #############################################
846
1037        def __stoppingcriterium_GMRES(self,norm_r,norm_b):        def __stoppingcriterium_GMRES(self,norm_r,norm_b):
1038            return self.stoppingcriterium_GMRES(norm_r,norm_b)            return self.stoppingcriterium_GMRES(norm_r,norm_b)
1039
1040          def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
1041              return self.stoppingcriterium_MINRES(norm_r,norm_Ax)
1042
1043
1044        def setTolerance(self,tolerance=1.e-8):        def setTolerance(self,tolerance=1.e-8):
1045                self.__tol=tolerance                self.__tol=tolerance
1046        def getTolerance(self):        def getTolerance(self):
1084                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1085                  #       u=v+(u-v)                  #       u=v+(u-v)
1086          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1087
1088              if solver=='MINRES':
1089                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1090                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
1091                    # solve Au=f-B^*p
1092                    #       A(u-v)=f-B^*p-Av
1093                    #       u=v+(u-v)
1094            u=v+self.solve_A(v,p)
1095
1096                else:                if solver=='PCG':
1097                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1098                  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)
1099              u=r[0]                u=r[0]
1100
1101                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1102
1103            return u,p            return u,p

Legend:
 Removed from v.1480 changed lines Added in v.1481