/[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 1518 by artak, Fri Apr 18 02:36:37 2008 UTC revision 1519 by artak, Tue Apr 22 03:45:36 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 585  def GMRES(b, Aprod, Msolve, bilinearform Line 587  def GMRES(b, Aprod, Msolve, bilinearform
587     while True:     while True:
588        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
589        x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)        x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
       iter+=iter_restart  
590        if stopped: break        if stopped: break
591          iter+=iter_restart    
592     return x     return x
593    
594  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):
# Line 613  def GMRESm(b, Aprod, Msolve, bilinearfor Line 615  def GMRESm(b, Aprod, Msolve, bilinearfor
615        
616     v.append(r/rho)     v.append(r/rho)
617     g[0]=rho     g[0]=rho
   
618     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     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      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
# Line 691  def GMRESm(b, Aprod, Msolve, bilinearfor Line 692  def GMRESm(b, Aprod, Msolve, bilinearfor
692    
693     return x,stopped     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):  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
831    
832      #      #
# Line 844  def MINRES(b, Aprod, Msolve, bilinearfor Line 980  def MINRES(b, Aprod, Msolve, bilinearfor
980    
981      return x      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):  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1050    
# Line 1139  class HomogeneousSaddlePointProblem(obje Line 1340  class HomogeneousSaddlePointProblem(obje
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'):        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1344            return self.stoppingcriterium2(norm_r,norm_b,solver)            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
# Line 1186  class HomogeneousSaddlePointProblem(obje Line 1387  class HomogeneousSaddlePointProblem(obje
1387                  #       u=v+(u-v)                  #       u=v+(u-v)
1388          u=v+self.solve_A(v,p)          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':                  if solver=='TFQMR':      
1400                  if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter                  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.)                  p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
# Line 1230  class HomogeneousSaddlePointProblem(obje Line 1440  class HomogeneousSaddlePointProblem(obje
1440        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1441            return self.B(v)            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     """     """
1451     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
Legend:
Removed from v.1518  
changed lines
  Added in v.1519
  ViewVC Help
Powered by ViewVC 1.1.26