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

revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1330 by gross, Mon Oct 22 04:54:49 2007 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} where C{norm_r}  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
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