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

revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1414 by gross, Thu Feb 14 10:01:43 2008 UTC
# Line 435  class NegativeNorm(SolverSchemeException Line 435  class NegativeNorm(SolverSchemeException
435     """     """
436     pass     pass
437
438  def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)):  class IterationHistory(object):
439       """
440       The IterationHistory class is used to define a stopping criterium. It keeps track of the
441       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
442       a given tolerance.
443       """
444       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
445          """
446          Initialization
447
448          @param tolerance: tolerance
449          @type tolerance: positive C{float}
450          @param verbose: switches on the printing out some information
451          @type verbose: C{bool}
452          """
453          if not tolerance>0.:
454              raise ValueError,"tolerance needs to be positive."
455          self.tolerance=tolerance
456          self.verbose=verbose
457          self.history=[]
458       def stoppingcriterium(self,norm_r,r,x):
459           """
460           returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.
461
462
463           @param norm_r: current residual norm
464           @type norm_r: non-negative C{float}
465           @param r: current residual (not used)
466           @param x: current solution approximation (not used)
467           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
468           @rtype: C{bool}
469
470           """
471           self.history.append(norm_r)
472           if self.verbose: print "iter: %s:  inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])
473           return self.history[-1]<=self.tolerance * self.history[0]
474
475    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
476     """     """
477     Solver for     Solver for
478
# Line 444  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 481  def PCG(b,x,Aprod,Msolve,bilinearform, n
481     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
482     It uses the conjugate gradient method with preconditioner M providing an approximation of A.     It uses the conjugate gradient method with preconditioner M providing an approximation of A.
483
484     The iteration is terminated if     The iteration is terminated if the C{stoppingcriterium} function return C{True}.

M{norm(r) <= tolerance * norm(b)}

where C{norm()} defines a norm and

M{r = b-Ax}

the residual.
485
486     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
487
# Line 461  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 490  def PCG(b,x,Aprod,Msolve,bilinearform, n
490     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst.
491
492     @param b: the right hand side of the liner system. C{b} is altered.     @param b: the right hand side of the liner system. C{b} is altered.
493     @type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y)     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
@param x: an initial guess for the solution
@type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y)
494     @param Aprod: returns the value Ax     @param Aprod: returns the value Ax
495     @type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R.     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.
496     @param Msolve: solves Mx=r     @param Msolve: solves Mx=r
497     @type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S.     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same
498    type like argument C{x}.
499     @param bilinearform: inner product C{<x,r>}     @param bilinearform: inner product C{<x,r>}
500     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}.     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.
501     @param norm: norm calculation for the residual C{r=b-Ax}.     @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.
502     @type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}.     @type stoppingcriterium: function that returns C{True} or C{False}
503     @param verbose: switches on the printing out some information     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
504     @type verbose: C{bool}     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
505     @param iter_max: maximum number of iteration steps.     @param iter_max: maximum number of iteration steps.
506     @type iter_max: C{int}     @type iter_max: C{int}
507     @param tolerance: tolerance     @return: the solution approximation and the corresponding residual
508     @type tolerance: positive C{float}     @rtype: C{tuple}
509     @return: the solution apprximation and the corresponding residual     @warning: C{b} and C{x} are altered.
@rtype: C{tuple} of an S type and and an R type object.A
@warning: C{b} ans C{x} are altered.
510     """     """
if verbose:
print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance   =%e"%(iter_max,tolerance)
511     iter=0     iter=0
512     normb = norm(b)     if x==None:
513     if normb<0: raise NegativeNorm        x=0*b
514       else:
515     b += (-1)*Aprod(x)        b += (-1)*Aprod(x)
516     r=b     r=b
517     rhat=Msolve(r)     rhat=Msolve(r)
518     d = rhat;     d = rhat
519     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
520       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521
522     while True:     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
523         normr=norm(r)         iter+=1
if normr<0: raise NegativeNorm
if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb)
if normr <= tolerance * normb: return x,r

iter+=1 # k = iter = 1 first time through
524         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
525
526         q=Aprod(d)         q=Aprod(d)
# Line 515  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 535  def PCG(b,x,Aprod,Msolve,bilinearform, n
535         d=rhat         d=rhat
536
537         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
538           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
539
540       return x,r
541
542    class ArithmeticTuple(object):
543       """
544       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
545
546       example of usage:
547
548       from esys.escript import Data
549       from numarray import array
550       a=Data(...)
551       b=array([1.,4.])
552       x=ArithmeticTuple(a,b)
553       y=5.*x
554
555       """
556       def __init__(self,*args):
557           """
558           initialize object with elements args.
559
560           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
561           """
562           self.__items=list(args)
563
564       def __len__(self):
565           """
566           number of items
567
568           @return: number of items
569           @rtype: C{int}
570           """
571           return len(self.__items)
572
573       def __getitem__(self,index):
574           """
575           get an item
576
577           @param index: item to be returned
578           @type index: C{int}
579           @return: item with index C{index}
580           """
581           return self.__items.__getitem__(index)
582
583       def __mul__(self,other):
584           """
585           scaling from the right
586
587           @param other: scaling factor
588           @type other: C{float}
589           @return: itemwise self*other
590           @rtype: L{ArithmeticTuple}
591           """
592           out=[]
593           for i in range(len(self)):
594               out.append(self[i]*other)
595           return ArithmeticTuple(*tuple(out))
596
597       def __rmul__(self,other):
598           """
599           scaling from the left
600
601           @param other: scaling factor
602           @type other: C{float}
603           @return: itemwise other*self
604           @rtype: L{ArithmeticTuple}
605           """
606           out=[]
607           for i in range(len(self)):
608               out.append(other*self[i])
609           return ArithmeticTuple(*tuple(out))
610
612           """
613           in-place add of other to self
614
615           @param other: increment
616           @type other: C{ArithmeticTuple}
617           """
618           if len(self) != len(other):
619               raise ValueError,"tuple length must match."
620           for i in range(len(self)):
621               self.__items[i]+=other[i]
622           return self
623
625          """
626          This provides a framwork for solving homogeneous saddle point problem of the form
627
628                 Av+B^*p=f
629                 Bv    =0
630
631          for the unknowns v and p and given operators A and B and given right hand side f.
632          B^* is the adjoint operator of B is the given inner product.
633
634          """
635          def __init__(self,**kwargs):
636            self.setTolerance()
637            self.setToleranceReductionFactor()
638
639          def initialize(self):
640            """
641            initialize the problem (overwrite)
642            """
643            pass
644          def B(self,v):
645             """
646             returns Bv (overwrite)
647             @rtype: equal to the type of p
648
649             @note: boundary conditions on p should be zero!
650             """
651             pass
652
653          def inner(self,p0,p1):
654             """
655             returns inner product of two element p0 and p1  (overwrite)
656
657             @type p0: equal to the type of p
658             @type p1: equal to the type of p
659             @rtype: C{float}
660
661             @rtype: equal to the type of p
662             """
663             pass
664
665          def solve_A(self,u,p):
666             """
667             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
668
669             @rtype: equal to the type of v
670             @note: boundary conditions on v should be zero!
671             """
672             pass
673
674          def solve_prec(self,p):
675             """
676             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
677
678             @rtype: equal to the type of p
679             """
680             pass
681
682          def stoppingcriterium(self,Bv,v,p):
683             """
684             returns a True if iteration is terminated. (overwrite)
685
686             @rtype: C{bool}
687             """
688             pass
689
690          def __inner(self,p,r):
691             return self.inner(p,r[1])
692
693          def __stoppingcriterium(self,norm_r,r,p):
694              return self.stoppingcriterium(r[1],r[0],p)
695
696          def setTolerance(self,tolerance=1.e-8):
697                  self.__tol=tolerance
698          def getTolerance(self):
699                  return self.__tol
700          def setToleranceReductionFactor(self,reduction=0.01):
701                  self.__reduction=reduction
702          def getSubProblemTolerance(self):
703                  return self.__reduction*self.getTolerance()
704
705          def solve(self,v,p,max_iter=20, verbose=False, show_details=False):
706                  """
707                  solves the saddle point problem using initial guesses v and p.
708
709                  @param max_iter: maximum number of iteration steps.
710                  """
711                  self.verbose=verbose
712                  self.show_details=show_details and self.verbose
713
714              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
715
716              self.__z=v+self.solve_A(v,p*0)
717                  Bz=self.B(self.__z)
718                  #
719              #   solve BA^-1B^*p = Bz
720                  #
721                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
722                  #
723                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
724                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
725                  #
726                  self.iter=0
727                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
728                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
729              return r[0],p
730
731          def __Msolve(self,r):
732              return self.solve_prec(r[1])
733
734          def __Aprod(self,p):
735              # return BA^-1B*p
736              #solve Av =-B^*p as Av =f-Az-B^*p
737              v=self.solve_A(self.__z,p)
738              return ArithmeticTuple(v, self.B(v))
739
740