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

revision 1474 by gross, Thu Apr 3 05:16:56 2008 UTC revision 1475 by artak, Mon Apr 7 05:13:21 2008 UTC
# Line 579  def givapp(c,s,vin): Line 579  def givapp(c,s,vin):
579              vrot[i:i+2]=w1,w2              vrot[i:i+2]=w1,w2
580      return vrot      return vrot
581
582  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
583       m=iter_restart
584       iter=0
585       while True:
586          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
587          x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
588          iter+=iter_restart
589          if stopped: break
590       return x
591
592    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
593     iter=0     iter=0
594     r=Msolve(b)     r=Msolve(b)
595     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
# Line 593  def GMRES(b, Aprod, Msolve, bilinearform Line 603  def GMRES(b, Aprod, Msolve, bilinearform
603        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
604        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm,"negative norm."
605
606     h=numarray.zeros((iter_max,iter_max),numarray.Float64)     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
607     c=numarray.zeros(iter_max,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
608     s=numarray.zeros(iter_max,numarray.Float64)     s=numarray.zeros(iter_restart,numarray.Float64)
609     g=numarray.zeros(iter_max,numarray.Float64)     g=numarray.zeros(iter_restart,numarray.Float64)
610     v=[]     v=[]
611
612     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
613     v.append(r/rho)     v.append(r/rho)
614     g[0]=rho     g[0]=rho
615
616     while not stoppingcriterium(rho,norm_b):     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
617
618      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
619
# Line 673  def GMRES(b, Aprod, Msolve, bilinearform Line 683  def GMRES(b, Aprod, Msolve, bilinearform
683     else : xhat=v[0]     else : xhat=v[0]
684
685     x += xhat     x += xhat
686       if iter!=iter_restart-1:
687          stopped=True
688       else:
689          stopped=False
690
691     return x     return x,stopped
692
693  #############################################  #############################################
694