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

revision 1554 by artak, Fri May 9 02:50:49 2008 UTC revision 1557 by artak, Mon May 19 04:44:27 2008 UTC
# Line 618  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
622     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
623
624      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
625

626      p=Msolve(Aprod(v[iter]))      p=Msolve(Aprod(v[iter]))
627
628      v.append(p)      v.append(p)
# Line 630  def GMRESm(b, Aprod, Msolve, bilinearfor Line 630  def GMRESm(b, Aprod, Msolve, bilinearfor
630      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
631
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]        v[iter+1]-=h[j][iter]*v[j]
636
637      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]))
638      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1][iter]
639

640  # Reorthogonalize if needed  # Reorthogonalize if needed
641      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
642       for j in range(iter+1):       for j in range(iter+1):
643          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
644              h[j][iter]=h[j][iter]+hr #vhat              h[j][iter]=h[j][iter]+hr
645              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
646
647       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
# Line 660  def GMRESm(b, Aprod, Msolve, bilinearfor Line 659  def GMRESm(b, Aprod, Msolve, bilinearfor
659              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i][iter]=hhat[i]
660
661      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
662
663      if mu!=0 :      if mu!=0 :
664          c[iter]=h[iter][iter]/mu          c[iter]=h[iter][iter]/mu
665          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1][iter]/mu
# Line 668  def GMRESm(b, Aprod, Msolve, bilinearfor Line 668  def GMRESm(b, Aprod, Msolve, bilinearfor
668          g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])          g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
669
670  # Update the residual norm  # Update the residual norm
671
672          rho=abs(g[iter+1])          rho=abs(g[iter+1])
673      iter+=1      iter+=1
674
# Line 686  def GMRESm(b, Aprod, Msolve, bilinearfor Line 687  def GMRESm(b, Aprod, Msolve, bilinearfor
687       for i in range(iter-1):       for i in range(iter-1):
688      xhat += v[i]*y[i]      xhat += v[i]*y[i]
689     else : xhat=v[0]     else : xhat=v[0]
690
691     x += xhat     x += xhat
692     if iter<iter_restart-1:     if iter<iter_restart-1:
693        stopped=True        stopped=True
# Line 1204  class ArithmeticTuple(object): Line 1205  class ArithmeticTuple(object):
1205         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1206         """         """
1207         out=[]         out=[]
1208         for i in range(len(self)):         other=1.*other
1209           if isinstance(other,float):
1210        for i in range(len(self)):
1211             out.append(self[i]*other)             out.append(self[i]*other)
1212           else:
1213            for i in range(len(self)):
1214               out.append(self[i]*other[i])
1215         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1216
1217     def __rmul__(self,other):     def __rmul__(self,other):
# Line 1218  class ArithmeticTuple(object): Line 1224  class ArithmeticTuple(object):
1224         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1225         """         """
1226         out=[]         out=[]
1227         for i in range(len(self)):         other=1.*other
1228           if isinstance(other,float):
1229        for i in range(len(self)):
1230             out.append(other*self[i])             out.append(other*self[i])
1231           else:
1232            for i in range(len(self)):
1233               out.append(other[i]*self[i])
1234         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1235
1236  #########################  #########################
# Line 1235  class ArithmeticTuple(object): Line 1246  class ArithmeticTuple(object):
1246         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1247         """         """
1248         out=[]         out=[]
1249         for i in range(len(self)):         other=1.*other
1250           if isinstance(other,float):
1251        for i in range(len(self)):
1252             out.append(self[i]*(1./other))             out.append(self[i]*(1./other))
1253           else:
1254            for i in range(len(self)):
1255               out.append(self[i]*(1./other[i]))
1256         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1257
1258     def __rdiv__(self,other):     def __rdiv__(self,other):
# Line 1249  class ArithmeticTuple(object): Line 1265  class ArithmeticTuple(object):
1265         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1266         """         """
1267         out=[]         out=[]
1268         for i in range(len(self)):         other=1.*other
1269             out.append(other/self[i])         if isinstance(other,float):
1270            for i in range(len(self)):
1271               out.append(other*(1./self[i]))
1272           else:
1273            for i in range(len(self)):
1274               out.append(other[i]*(1./self[i]))
1275         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1276
1277  ##########################################33  ##########################################33
# Line 1298  class ArithmeticTuple(object): Line 1319  class ArithmeticTuple(object):
1319         for i in range(len(self)):         for i in range(len(self)):
1320             self.__items[i]-=other[i]             self.__items[i]-=other[i]
1321         return self         return self
1322
1323       def __isub__(self,other):
1324           """
1325           subtract other from self
1326
1327           @param other: increment
1328           @type other: C{ArithmeticTuple}
1329           """
1330           if len(self) != len(other):
1331               raise ValueError,"tuple length must match."
1332           for i in range(len(self)):
1333               self.__items[i]-=other[i]
1334           return self
1335
1336     def __neg__(self):     def __neg__(self):
1337         """         """
# Line 1416  class HomogeneousSaddlePointProblem(obje Line 1450  class HomogeneousSaddlePointProblem(obje
1450
1451            # 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)
1452            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
1453                Bz=self.B(self.__z)                Bz=self.B(self.__z)
1454                #                #
1455            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1456                #                #
# Line 1458  class HomogeneousSaddlePointProblem(obje Line 1492  class HomogeneousSaddlePointProblem(obje
1492
1493            if solver=='GMRESC':                  if solver=='GMRESC':
1494                  if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter                  if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1495                  p0=self.solve_prec(Bz)                  p0=self.solve_prec1(Bz)
1496             #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))              #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1497                 #p-=alfa                  #p-=alfa
1498                  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)                  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)
1499                  #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())                  #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())
1500
1501                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1502                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1503                  #       u=v+(u-v)                  #       u=v+(u-v)
1504              p=x[1]              p=x[1]
1505          #u=v+self.solve_A(v,p)                u=v+self.solve_A(v,p)
1506          #p=x[1]          #p=x[1]
1507          #u=x[0]          #u=x[0]
1508

Legend:
 Removed from v.1554 changed lines Added in v.1557

 ViewVC Help Powered by ViewVC 1.1.26