/[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 1956 by gross, Mon Nov 3 05:08:42 2008 UTC revision 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Provides some tools related to PDEs.  Provides some tools related to PDEs.
24    
25  Currently includes:  Currently includes:
26      - Projector - to project a discontinuous      - Projector - to project a discontinuous function onto a continuous function
27      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
28      - TimeIntegrationManager - to handel extraplotion in time      - TimeIntegrationManager - to handle extrapolation in time
29      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30    
31  @var __author__: name of author  :var __author__: name of author
32  @var __copyright__: copyrights  :var __copyright__: copyrights
33  @var __license__: licence agreement  :var __license__: licence agreement
34  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
35  @var __version__: version  :var __version__: version
36  @var __date__: date of the version  :var __date__: date of the version
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
# Line 41  __author__="Lutz Gross, l.gross@uq.edu.a Line 41  __author__="Lutz Gross, l.gross@uq.edu.a
41    
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numpy
45  import util  import util
46  import math  import math
47    
 ##### Added by Artak  
 # from Numeric import zeros,Int,Float64  
 ###################################  
   
   
48  class TimeIntegrationManager:  class TimeIntegrationManager:
49    """    """
50    a simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
51    
52    typical usage is::    Typical usage is::
53    
54       dt=0.1 # time increment       dt=0.1 # time increment
55       tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
# Line 64  class TimeIntegrationManager: Line 59  class TimeIntegrationManager:
59           tm.checkin(dt,v)           tm.checkin(dt,v)
60           t+=dt           t+=dt
61    
62    @note: currently only p=1 is supported.    :note: currently only p=1 is supported.
63    """    """
64    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
65       """       """
66       sets up the value manager where inital_value is the initial value and p is order used for extrapolation       Sets up the value manager where ``inital_values`` are the initial values
67         and p is the order used for extrapolation.
68       """       """
69       if kwargs.has_key("p"):       if kwargs.has_key("p"):
70              self.__p=kwargs["p"]              self.__p=kwargs["p"]
# Line 85  class TimeIntegrationManager: Line 81  class TimeIntegrationManager:
81    
82    def getTime(self):    def getTime(self):
83        return self.__t        return self.__t
84    
85    def getValue(self):    def getValue(self):
86        out=self.__v_mem[0]        out=self.__v_mem[0]
87        if len(out)==1:        if len(out)==1:
# Line 94  class TimeIntegrationManager: Line 91  class TimeIntegrationManager:
91    
92    def checkin(self,dt,*values):    def checkin(self,dt,*values):
93        """        """
94        adds new values to the manager. the p+1 last value get lost        Adds new values to the manager. The p+1 last values are lost.
95        """        """
96        o=min(self.__order+1,self.__p)        o=min(self.__order+1,self.__p)
97        self.__order=min(self.__order+1,self.__p)        self.__order=min(self.__order+1,self.__p)
# Line 111  class TimeIntegrationManager: Line 108  class TimeIntegrationManager:
108    
109    def extrapolate(self,dt):    def extrapolate(self,dt):
110        """        """
111        extrapolates to dt forward in time.        Extrapolates to ``dt`` forward in time.
112        """        """
113        if self.__order==0:        if self.__order==0:
114           out=self.__v_mem[0]           out=self.__v_mem[0]
# Line 126  class TimeIntegrationManager: Line 123  class TimeIntegrationManager:
123           return out[0]           return out[0]
124        else:        else:
125           return out           return out
126    
127    
128  class Projector:  class Projector:
129    """    """
130    The Projector is a factory which projects a discontiuous function onto a    The Projector is a factory which projects a discontinuous function onto a
131    continuous function on the a given domain.    continuous function on a given domain.
132    """    """
133    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce=True, fast=True):
134      """      """
135      Create a continuous function space projector for a domain.      Creates a continuous function space projector for a domain.
136    
137      @param domain: Domain of the projection.      :param domain: Domain of the projection.
138      @param reduce: Flag to reduce projection order (default is True)      :param reduce: Flag to reduce projection order
139      @param fast: Flag to use a fast method based on matrix lumping (default is true)      :param fast: Flag to use a fast method based on matrix lumping
140      """      """
141      self.__pde = linearPDEs.LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
142      if fast:      if fast:
143        self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)          self.__pde.getSolverOptions().setSolverMethod(linearPDEs.SolverOptions.LUMPING)
144      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
145      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
146      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
147      return      return
148      def getSolverOptions(self):
149        """
150        Returns the solver options of the PDE solver.
151        
152        :rtype: `linearPDEs.SolverOptions`
153        """
154        return self.__pde.getSolverOptions()
155    
156      def getValue(self, input_data):
157        """
158        Projects ``input_data`` onto a continuous function.
159    
160        :param input_data: the data to be projected
161        """
162        return self(input_data)
163    
164    def __call__(self, input_data):    def __call__(self, input_data):
165      """      """
166      Projects input_data onto a continuous function      Projects ``input_data`` onto a continuous function.
167    
168      @param input_data: The input_data to be projected.      :param input_data: the data to be projected
169      """      """
170      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
171      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
# Line 186  class Projector: Line 198  class Projector:
198    
199  class NoPDE:  class NoPDE:
200       """       """
201       solves the following problem for u:       Solves the following problem for u:
202    
203       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
204    
205       with constraint       with constraint
206    
207       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
208    
209       where D, Y, r and q are given functions of rank 1.       where *D*, *Y*, *r* and *q* are given functions of rank 1.
210    
211       In the case of scalars this takes the form       In the case of scalars this takes the form
212    
213       M{D*u=Y}       *D*u=Y*
214    
215       with constraint       with constraint
216    
217       M{u=r}  where M{q>0}       *u=r* where *q>0*
218    
219       where D, Y, r and q are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
220    
221       The constraint is overwriting any other condition.       The constraint overwrites any other condition.
222    
223       @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention       :note: This class is similar to the `linearPDEs.LinearPDE` class with
224              that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole              A=B=C=X=0 but has the intention that all input parameters are given
225              thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.              in `Solution` or `ReducedSolution`.
226       """       """
227         # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
228         # this.
229       def __init__(self,domain,D=None,Y=None,q=None,r=None):       def __init__(self,domain,D=None,Y=None,q=None,r=None):
230           """           """
231           initialize the problem           Initializes the problem.
232    
233           @param domain: domain of the PDE.           :param domain: domain of the PDE
234           @type domain: L{Domain}           :type domain: `Domain`
235           @param D: coefficient of the solution.           :param D: coefficient of the solution
236           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
237           @param Y: right hand side           :param Y: right hand side
238           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
239           @param q: location of constraints           :param q: location of constraints
240           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
241           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
242           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
243           """           """
244           self.__domain=domain           self.__domain=domain
245           self.__D=D           self.__D=D
# Line 234  class NoPDE: Line 248  class NoPDE:
248           self.__r=r           self.__r=r
249           self.__u=None           self.__u=None
250           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
251    
252       def setReducedOn(self):       def setReducedOn(self):
253           """           """
254           sets the L{FunctionSpace} of the solution to L{ReducedSolution}           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
255           """           """
256           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
257           self.__u=None           self.__u=None
258    
259       def setReducedOff(self):       def setReducedOff(self):
260           """           """
261           sets the L{FunctionSpace} of the solution to L{Solution}           Sets the `FunctionSpace` of the solution to `Solution`.
262           """           """
263           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
264           self.__u=None           self.__u=None
265            
266       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
267           """           """
268           assigns values to the parameters.           Assigns values to the parameters.
269    
270           @param D: coefficient of the solution.           :param D: coefficient of the solution
271           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
272           @param Y: right hand side           :param Y: right hand side
273           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
274           @param q: location of constraints           :param q: location of constraints
275           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
276           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
277           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
278           """           """
279           if not D==None:           if not D==None:
280              self.__D=D              self.__D=D
# Line 276  class NoPDE: Line 291  class NoPDE:
291    
292       def getSolution(self):       def getSolution(self):
293           """           """
294           returns the solution           Returns the solution.
295            
296           @return: the solution of the problem           :return: the solution of the problem
297           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.           :rtype: `Data` object in the `FunctionSpace` `Solution` or
298                     `ReducedSolution`
299           """           """
300           if self.__u==None:           if self.__u==None:
301              if self.__D==None:              if self.__D==None:
# Line 296  class NoPDE: Line 312  class NoPDE:
312                  self.__u*=(1.-q)                  self.__u*=(1.-q)
313                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
314           return self.__u           return self.__u
315                
316  class Locator:  class Locator:
317       """       """
318       Locator provides access to the values of data objects at a given       Locator provides access to the values of data objects at a given spatial
319       spatial coordinate x.         coordinate x.
320        
321       In fact, a Locator object finds the sample in the set of samples of a       In fact, a Locator object finds the sample in the set of samples of a
322       given function space or domain where which is closest to the given       given function space or domain which is closest to the given point x.
      point x.  
323       """       """
324    
325       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
326         """         """
327         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
328         or FunctionSpace where for the sample point which         or FunctionSpace for the sample point which is closest to the given
329         closest to the given point x.         point x.
330    
331         @param where: function space         :param where: function space
332         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
333         @param x: coefficient of the solution.         :param x: location(s) of the Locator
334         @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
335         """         """
336         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
337            self.__function_space=where            self.__function_space=where
338         else:         else:
339            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
340           iterative=False
341         if isinstance(x, list):         if isinstance(x, list):
342               if len(x)==0:
343                  raise "ValueError", "At least one point must be given."
344               try:
345                 iter(x[0])
346                 iterative=True
347               except TypeError:
348                 iterative=False
349           if iterative:
350             self.__id=[]             self.__id=[]
351             for p in x:             for p in x:
352                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
# Line 334  class Locator: Line 358  class Locator:
358         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
359         """         """
360         x=self.getX()         x=self.getX()
361         if instance(x,list):         if isinstance(x,list):
362            out="["            out="["
363            first=True            first=True
364            for xx in x:            for xx in x:
# Line 350  class Locator: Line 374  class Locator:
374    
375       def getX(self):       def getX(self):
376          """          """
377      Returns the exact coordinates of the Locator.          Returns the exact coordinates of the Locator.
378      """          """
379          return self(self.getFunctionSpace().getX())          return self(self.getFunctionSpace().getX())
380    
381       def getFunctionSpace(self):       def getFunctionSpace(self):
382          """          """
383      Returns the function space of the Locator.          Returns the function space of the Locator.
384      """          """
385          return self.__function_space          return self.__function_space
386    
387       def getId(self,item=None):       def getId(self,item=None):
388          """          """
389      Returns the identifier of the location.          Returns the identifier of the location.
390      """          """
391          if item == None:          if item == None:
392             return self.__id             return self.__id
393          else:          else:
# Line 375  class Locator: Line 399  class Locator:
399    
400       def __call__(self,data):       def __call__(self,data):
401          """          """
402      Returns the value of data at the Locator of a Data object otherwise          Returns the value of data at the Locator of a Data object.
403      the object is returned.          """
     """  
404          return self.getValue(data)          return self.getValue(data)
405    
406       def getValue(self,data):       def getValue(self,data):
407          """          """
408      Returns the value of data at the Locator if data is a Data object          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
409      otherwise the object is returned.          object otherwise the object is returned.
410      """          """
411          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
412             if data.getFunctionSpace()==self.getFunctionSpace():             dat=util.interpolate(data,self.getFunctionSpace())
              dat=data  
            else:  
              dat=data.interpolate(self.getFunctionSpace())  
413             id=self.getId()             id=self.getId()
414             r=data.getRank()             r=data.getRank()
415             if isinstance(id,list):             if isinstance(id,list):
416                 out=[]                 out=[]
417                 for i in id:                 for i in id:
418                    o=data.getValueOfGlobalDataPoint(*i)                    o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
419                    if data.getRank()==0:                    if data.getRank()==0:
420                       out.append(o[0])                       out.append(o[0])
421                    else:                    else:
422                       out.append(o)                       out.append(o)
423                 return out                 return out
424             else:             else:
425               out=data.getValueOfGlobalDataPoint(*id)               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
426               if data.getRank()==0:               if data.getRank()==0:
427                  return out[0]                  return out[0]
428               else:               else:
# Line 410  class Locator: Line 430  class Locator:
430          else:          else:
431             return data             return data
432    
433    
434    def getInfLocator(arg):
435        """
436        Return a Locator for a point with the inf value over all arg.
437        """
438        if not isinstance(arg, escript.Data):
439        raise TypeError,"getInfLocator: Unknown argument type."
440        a_inf=util.inf(arg)
441        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
442        x=arg.getFunctionSpace().getX()
443        x_min=x.getTupleForGlobalDataPoint(*loc)
444        return Locator(arg.getFunctionSpace(),x_min)
445    
446    def getSupLocator(arg):
447        """
448        Return a Locator for a point with the sup value over all arg.
449        """
450        if not isinstance(arg, escript.Data):
451        raise TypeError,"getInfLocator: Unknown argument type."
452        a_inf=util.sup(arg)
453        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
454        x=arg.getFunctionSpace().getX()
455        x_min=x.getTupleForGlobalDataPoint(*loc)
456        return Locator(arg.getFunctionSpace(),x_min)
457        
458    
459  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
460     """     """
461     exceptions thrown by solvers     This is a generic exception thrown by solvers.
462     """     """
463     pass     pass
464    
465  class IndefinitePreconditioner(SolverSchemeException):  class IndefinitePreconditioner(SolverSchemeException):
466     """     """
467     the preconditioner is not positive definite.     Exception thrown if the preconditioner is not positive definite.
468     """     """
469     pass     pass
470    
471  class MaxIterReached(SolverSchemeException):  class MaxIterReached(SolverSchemeException):
472     """     """
473     maxium number of iteration steps is reached.     Exception thrown if the maximum number of iteration steps is reached.
474     """     """
475     pass     pass
476  class IterationBreakDown(SolverSchemeException):  
477    class CorrectionFailed(SolverSchemeException):
478     """     """
479     iteration scheme econouters an incurable breakdown.     Exception thrown if no convergence has been achieved in the solution
480       correction scheme.
481     """     """
482     pass     pass
483  class NegativeNorm(SolverSchemeException):  
484    class IterationBreakDown(SolverSchemeException):
485     """     """
486     a norm calculation returns a negative norm.     Exception thrown if the iteration scheme encountered an incurable breakdown.
487     """     """
488     pass     pass
489    
490  class IterationHistory(object):  class NegativeNorm(SolverSchemeException):
491     """     """
492     The IterationHistory class is used to define a stopping criterium. It keeps track of the     Exception thrown if a norm calculation returns a negative norm.
    residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by  
    a given tolerance.  
493     """     """
494     def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):     pass
       """  
       Initialization  
   
       @param tolerance: tolerance  
       @type tolerance: positive C{float}  
       @param verbose: switches on the printing out some information  
       @type verbose: C{bool}  
       """  
       if not tolerance>0.:  
           raise ValueError,"tolerance needs to be positive."  
       self.tolerance=tolerance  
       self.verbose=verbose  
       self.history=[]  
    def stoppingcriterium(self,norm_r,r,x):  
        """  
        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.  
   
         
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param r: current residual (not used)  
        @param x: current solution approximation (not used)  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
   
        """  
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=self.tolerance * self.history[0]  
495    
496     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):  def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
497         """     """
498         returns True if the C{norm_r} is C{tolerance}*C{norm_b}     Solver for
499    
500             *Ax=b*
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param norm_b: norm of right hand side  
        @type norm_b: non-negative C{float}  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
501    
502         """     with a symmetric and positive definite operator A (more details required!).
503         if TOL==None:     It uses the conjugate gradient method with preconditioner M providing an
504            TOL=self.tolerance     approximation of A.
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=TOL * norm_b  
505    
506  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):     The iteration is terminated if
    """  
    Solver for  
507    
508     M{Ax=b}     *|r| <= atol+rtol*|r0|*
509    
510     with a symmetric and positive definite operator A (more details required!).     where *r0* is the initial residual and *|.|* is the energy norm. In fact
    It uses the conjugate gradient method with preconditioner M providing an approximation of A.  
511    
512     The iteration is terminated if the C{stoppingcriterium} function return C{True}.     *|r| = sqrt( bilinearform(Msolve(r),r))*
513    
514     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
515    
516     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,     I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
517     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
518     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst}.
519    
520     @param b: the right hand side of the liner system. C{b} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
521     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
522     @param Aprod: returns the value Ax     :param x: an initial guess for the solution
523     @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}.     :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
524     @param Msolve: solves Mx=r     :param Aprod: returns the value Ax
525     @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     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
526  type like argument C{x}.                  argument ``x``. The returned object needs to be of the same type
527     @param bilinearform: inner product C{<x,r>}                  like argument ``r``.
528     @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}.     :param Msolve: solves Mx=r
529     @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.     :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
530     @type stoppingcriterium: function that returns C{True} or C{False}                   argument ``r``. The returned object needs to be of the same
531     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.                   type like argument ``x``.
532     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :param bilinearform: inner product ``<x,r>``
533     @param iter_max: maximum number of iteration steps.     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
534     @type iter_max: C{int}                         type like argument ``x`` and ``r`` is. The returned value
535     @return: the solution approximation and the corresponding residual                         is a ``float``.
536     @rtype: C{tuple}     :param atol: absolute tolerance
537     @warning: C{b} and C{x} are altered.     :type atol: non-negative ``float``
538       :param rtol: relative tolerance
539       :type rtol: non-negative ``float``
540       :param iter_max: maximum number of iteration steps
541       :type iter_max: ``int``
542       :return: the solution approximation and the corresponding residual
543       :rtype: ``tuple``
544       :warning: ``r`` and ``x`` are altered.
545     """     """
546     iter=0     iter=0
    if x==None:  
       x=0*b  
    else:  
       b += (-1)*Aprod(x)  
    r=b  
547     rhat=Msolve(r)     rhat=Msolve(r)
548     d = rhat     d = rhat
549     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
550     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm,"negative norm."
551       norm_r0=math.sqrt(rhat_dot_r)
552       atol2=atol+rtol*norm_r0
553       if atol2<=0:
554          raise ValueError,"Non-positive tolarance."
555       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
556    
557       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
558    
559    
560     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     while not math.sqrt(rhat_dot_r) <= atol2:
561         iter+=1         iter+=1
562         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
563    
564         q=Aprod(d)         q=Aprod(d)
565         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
566         x += alpha * d         x += alpha * d
567         r += (-alpha) * q         if isinstance(q,ArithmeticTuple):
568           r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
569           else:
570               r += (-alpha) * q
571         rhat=Msolve(r)         rhat=Msolve(r)
572         rhat_dot_r_new = bilinearform(rhat, r)         rhat_dot_r_new = bilinearform(rhat, r)
573         beta = rhat_dot_r_new / rhat_dot_r         beta = rhat_dot_r_new / rhat_dot_r
# Line 557  type like argument C{x}. Line 576  type like argument C{x}.
576    
577         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
578         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
579           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
580     return x,r     if verbose: print "PCG: tolerance reached after %s steps."%iter
581       return x,r,math.sqrt(rhat_dot_r)
582    
583  class Defect(object):  class Defect(object):
584      """      """
585      defines a non-linear defect F(x) of a variable x      Defines a non-linear defect F(x) of a variable x.
586      """      """
587      def __init__(self):      def __init__(self):
588          """          """
589          initialize defect          Initializes defect.
590          """          """
591          self.setDerivativeIncrementLength()          self.setDerivativeIncrementLength()
592    
593      def bilinearform(self, x0, x1):      def bilinearform(self, x0, x1):
594          """          """
595          returns the inner product of x0 and x1          Returns the inner product of x0 and x1
596          @param x0: a value for x  
597          @param x1: a value for x          :param x0: value for x0
598          @return: the inner product of x0 and x1          :param x1: value for x1
599          @rtype: C{float}          :return: the inner product of x0 and x1
600            :rtype: ``float``
601          """          """
602          return 0          return 0
603          
604      def norm(self,x):      def norm(self,x):
605          """          """
606          the norm of argument C{x}          Returns the norm of argument ``x``.
607    
608          @param x: a value for x          :param x: a value
609          @return: norm of argument x          :return: norm of argument x
610          @rtype: C{float}          :rtype: ``float``
611          @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
612          """          """
613          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
614          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
615          return math.sqrt(s)          return math.sqrt(s)
616    
   
617      def eval(self,x):      def eval(self,x):
618          """          """
619          returns the value F of a given x          Returns the value F of a given ``x``.
620    
621          @param x: value for which the defect C{F} is evalulated.          :param x: value for which the defect ``F`` is evaluated
622          @return: value of the defect at C{x}          :return: value of the defect at ``x``
623          """          """
624          return 0          return 0
625    
626      def __call__(self,x):      def __call__(self,x):
627          return self.eval(x)          return self.eval(x)
628    
629      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
630          """          """
631          sets the relative length of the increment used to approximate the derivative of the defect          Sets the relative length of the increment used to approximate the
632          the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
633            direction of v with x as a starting point.
634    
635          @param inc: relative increment length          :param inc: relative increment length
636          @type inc: positive C{float}          :type inc: positive ``float``
637          """          """
638          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError,"positive increment required."
639          self.__inc=inc          self.__inc=inc
640    
641      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
642          """          """
643          returns the relative increment length used to approximate the derivative of the defect          Returns the relative increment length used to approximate the
644          @return: value of the defect at C{x}          derivative of the defect.
645          @rtype: positive C{float}          :return: value of the defect at ``x``
646            :rtype: positive ``float``
647          """          """
648          return self.__inc          return self.__inc
649    
650      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
651          """          """
652          returns the directional derivative at x0 in the direction of v          Returns the directional derivative at ``x0`` in the direction of ``v``.
653    
654          @param F0: value of this defect at x0          :param F0: value of this defect at x0
655          @param x0: value at which derivative is calculated.          :param x0: value at which derivative is calculated
656          @param v: direction          :param v: direction
657          @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)          :param v_is_normalised: True to indicate that ``v`` is nomalized
658          @return: derivative of this defect at x0 in the direction of C{v}                                  (self.norm(v)=0)
659          @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method          :return: derivative of this defect at x0 in the direction of ``v``
660          maybe oepsnew verwritten to use exact evalution.          :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
661                   used but this method maybe overwritten to use exact evaluation.
662          """          """
663          normx=self.norm(x0)          normx=self.norm(x0)
664          if normx>0:          if normx>0:
# Line 651  class Defect(object): Line 674  class Defect(object):
674          F1=self.eval(x0 + epsnew * v)          F1=self.eval(x0 + epsnew * v)
675          return (F1-F0)/epsnew          return (F1-F0)/epsnew
676    
677  ######################################      ######################################
678  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, subtol_max=0.5, gamma=0.9, verbose=False):
679     """     """
680     solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:     Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
681       criterion:
682    
683     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
684      
685     where M{x0} is the initial guess.     where *x0* is the initial guess.
686    
687     @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.     :param defect: object defining the function *F*. ``defect.norm`` defines the
688     @type defect: L{Defect}                    *norm* used in the stopping criterion.
689     @param x: initial guess for the solution, C{x} is altered.     :type defect: `Defect`
690     @type x: any object type allowing basic operations such as  L{numarray.NumArray}, L{Data}     :param x: initial guess for the solution, ``x`` is altered.
691     @param iter_max: maximum number of iteration steps     :type x: any object type allowing basic operations such as
692     @type iter_max: positive C{int}              ``numpy.ndarray``, `Data`
693     @param sub_iter_max:     :param iter_max: maximum number of iteration steps
694     @type sub_iter_max:     :type iter_max: positive ``int``
695     @param atol: absolute tolerance for the solution     :param sub_iter_max: maximum number of inner iteration steps
696     @type atol: positive C{float}     :type sub_iter_max: positive ``int``
697     @param rtol: relative tolerance for the solution     :param atol: absolute tolerance for the solution
698     @type rtol: positive C{float}     :type atol: positive ``float``
699     @param gamma: tolerance safety factor for inner iteration     :param rtol: relative tolerance for the solution
700     @type gamma: positive C{float}, less than 1     :type rtol: positive ``float``
701     @param sub_tol_max: upper bound for inner tolerance.     :param gamma: tolerance safety factor for inner iteration
702     @type sub_tol_max: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
703     @return: an approximation of the solution with the desired accuracy     :param subtol_max: upper bound for inner tolerance
704     @rtype: same type as the initial guess.     :type subtol_max: positive ``float``, less than 1
705       :return: an approximation of the solution with the desired accuracy
706       :rtype: same type as the initial guess
707     """     """
708     lmaxit=iter_max     lmaxit=iter_max
709     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError,"atol needs to be non-negative."
710     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError,"rtol needs to be non-negative."
711     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
712     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma
713     if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max     if subtol_max<=0 or subtol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max
714    
715     F=defect(x)     F=defect(x)
716     fnrm=defect.norm(F)     fnrm=defect.norm(F)
717     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
718     sub_tol=sub_tol_max     subtol=subtol_max
719     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
720     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print "             tolerance = %e."%subtol
721     iter=1     iter=1
722     #     #
723     # main iteration loop     # main iteration loop
724     #     #
725     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
726              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
727              #              #
728          #   adjust sub_tol_          #   adjust subtol_
729          #          #
730              if iter > 1:              if iter > 1:
731             rat=fnrm/fnrmo             rat=fnrm/fnrmo
732                 sub_tol_old=sub_tol                 subtol_old=subtol
733             sub_tol=gamma*rat**2             subtol=gamma*rat**2
734             if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)             if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
735             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
736          #          #
737          # calculate newton increment xc          # calculate newton increment xc
738              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
739              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
740              #     if  atol is reached sub_iter returns the numer of steps performed to get there              #     if  atol is reached sub_iter returns the numer of steps performed to get there
741              #              #
742              #                #
743              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%subtol
744              try:              try:
745                 xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)                 xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
746              except MaxIterReached:              except MaxIterReached:
747                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
748              if sub_iter<0:              if sub_iter<0:
# Line 734  def NewtonGMRES(defect, x, iter_max=100, Line 760  def NewtonGMRES(defect, x, iter_max=100,
760    
761  def __givapp(c,s,vin):  def __givapp(c,s,vin):
762      """      """
763      apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin      Applies a sequence of Givens rotations (c,s) recursively to the vector
764      @warning: C{vin} is altered.      ``vin``
765    
766        :warning: ``vin`` is altered.
767      """      """
768      vrot=vin      vrot=vin
769      if isinstance(c,float):      if isinstance(c,float):
770          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
771      else:      else:
772          for i in range(len(c)):          for i in range(len(c)):
773              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
774          w2=s[i]*vrot[i]+c[i]*vrot[i+1]          w2=s[i]*vrot[i]+c[i]*vrot[i+1]
775              vrot[i:i+2]=w1,w2              vrot[i]=w1
776                vrot[i+1]=w2
777      return vrot      return vrot
778    
779  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
780     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
781     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
782     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
783     g=numarray.zeros(iter_restart,numarray.Float64)     g=numpy.zeros(iter_restart,numpy.float64)
784     v=[]     v=[]
785    
786     rho=defect.norm(F0)     rho=defect.norm(F0)
787     if rho<=0.: return x0*0     if rho<=0.: return x0*0
788      
789     v.append(-F0/rho)     v.append(-F0/rho)
790     g[0]=rho     g[0]=rho
791     iter=0     iter=0
792     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
793            if iter  >= iter_max:
794      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached,"maximum number of %s steps reached."%iter_max
795    
796          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
797      v.append(p)          v.append(p)
798    
799      v_norm1=defect.norm(v[iter+1])          v_norm1=defect.norm(v[iter+1])
800    
801          # Modified Gram-Schmidt          # Modified Gram-Schmidt
802      for j in range(iter+1):          for j in range(iter+1):
803           h[j][iter]=defect.bilinearform(v[j],v[iter+1])                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
804           v[iter+1]-=h[j][iter]*v[j]              v[iter+1]-=h[j,iter]*v[j]
805          
806      h[iter+1][iter]=defect.norm(v[iter+1])          h[iter+1,iter]=defect.norm(v[iter+1])
807      v_norm2=h[iter+1][iter]          v_norm2=h[iter+1,iter]
808    
809          # Reorthogonalize if needed          # Reorthogonalize if needed
810      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)
811          for j in range(iter+1):                for j in range(iter+1):
812             hr=defect.bilinearform(v[j],v[iter+1])                  hr=defect.bilinearform(v[j],v[iter+1])
813                 h[j][iter]=h[j][iter]+hr                  h[j,iter]=h[j,iter]+hr
814                 v[iter+1] -= hr*v[j]                  v[iter+1] -= hr*v[j]
815    
816          v_norm2=defect.norm(v[iter+1])              v_norm2=defect.norm(v[iter+1])
817          h[iter+1][iter]=v_norm2              h[iter+1,iter]=v_norm2
818          #   watch out for happy breakdown          #   watch out for happy breakdown
819          if not v_norm2 == 0:          if not v_norm2 == 0:
820                  v[iter+1]=v[iter+1]/h[iter+1][iter]              v[iter+1]=v[iter+1]/h[iter+1,iter]
821    
822          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
823      if iter > 0 :          if iter > 0 :
824          hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numpy.zeros(iter+1,numpy.float64)
825          for i in range(iter+1) : hhat[i]=h[i][iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
826          hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
827              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
828    
829      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])
830    
831      if mu!=0 :          if mu!=0 :
832          c[iter]=h[iter][iter]/mu              c[iter]=h[iter,iter]/mu
833          s[iter]=-h[iter+1][iter]/mu              s[iter]=-h[iter+1,iter]/mu
834          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]              h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
835          h[iter+1][iter]=0.0              h[iter+1,iter]=0.0
836          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])              gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
837                g[iter]=gg[0]
838                g[iter+1]=gg[1]
839    
840          # Update the residual norm          # Update the residual norm
841          rho=abs(g[iter+1])          rho=abs(g[iter+1])
842      iter+=1          iter+=1
843    
844     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
845     # It's time to compute x and leave.             # It's time to compute x and leave.
846     if iter > 0 :     if iter > 0 :
847       y=numarray.zeros(iter,numarray.Float64)           y=numpy.zeros(iter,numpy.float64)
848       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
849       if iter > 1 :         if iter > 1 :
850          i=iter-2            i=iter-2
851          while i>=0 :          while i>=0 :
852            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
853            i=i-1            i=i-1
854       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
855       for i in range(iter-1):       for i in range(iter-1):
856      xhat += v[i]*y[i]      xhat += v[i]*y[i]
857     else :     else :
858        xhat=v[0] * 0        xhat=v[0] * 0
859    
860     if iter<iter_restart-1:     if iter<iter_restart-1:
861        stopped=iter        stopped=iter
862     else:     else:
863        stopped=-1        stopped=-1
864    
865     return xhat,stopped     return xhat,stopped
866    
867  ##############################################  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
868  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):     """
869  ################################################     Solver for
870    
871       *Ax=b*
872    
873       with a general operator A (more details required!).
874       It uses the generalized minimum residual method (GMRES).
875    
876       The iteration is terminated if
877    
878       *|r| <= atol+rtol*|r0|*
879    
880       where *r0* is the initial residual and *|.|* is the energy norm. In fact
881    
882       *|r| = sqrt( bilinearform(r,r))*
883    
884       :param r: initial residual *r=b-Ax*. ``r`` is altered.
885       :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
886       :param x: an initial guess for the solution
887       :type x: same like ``r``
888       :param Aprod: returns the value Ax
889       :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
890                    argument ``x``. The returned object needs to be of the same
891                    type like argument ``r``.
892       :param bilinearform: inner product ``<x,r>``
893       :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
894                           type like argument ``x`` and ``r``. The returned value is
895                           a ``float``.
896       :param atol: absolute tolerance
897       :type atol: non-negative ``float``
898       :param rtol: relative tolerance
899       :type rtol: non-negative ``float``
900       :param iter_max: maximum number of iteration steps
901       :type iter_max: ``int``
902       :param iter_restart: in order to save memory the orthogonalization process
903                            is terminated after ``iter_restart`` steps and the
904                            iteration is restarted.
905       :type iter_restart: ``int``
906       :return: the solution approximation and the corresponding residual
907       :rtype: ``tuple``
908       :warning: ``r`` and ``x`` are altered.
909       """
910     m=iter_restart     m=iter_restart
911       restarted=False
912     iter=0     iter=0
913     xc=x     if rtol>0:
914          r_dot_r = bilinearform(r, r)
915          if r_dot_r<0: raise NegativeNorm,"negative norm."
916          atol2=atol+rtol*math.sqrt(r_dot_r)
917          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
918       else:
919          atol2=atol
920          if verbose: print "GMRES: absolute tolerance = %e"%atol2
921       if atol2<=0:
922          raise ValueError,"Non-positive tolarance."
923    
924     while True:     while True:
925        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
926        xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)        if restarted:
927             r2 = r-Aprod(x-x2)
928          else:
929             r2=1*r
930          x2=x*1.
931          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
932          iter+=iter_restart
933        if stopped: break        if stopped: break
934        iter+=iter_restart            if verbose: print "GMRES: restart."
935     return xc        restarted=True
936       if verbose: print "GMRES: tolerance has been reached."
937       return x
938    
939  def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
940     iter=0     iter=0
    r=Msolve(b)  
    r_dot_r = bilinearform(r, r)  
    if r_dot_r<0: raise NegativeNorm,"negative norm."  
    norm_b=math.sqrt(r_dot_r)  
941    
942     if x==None:     h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
943        x=0*b     c=numpy.zeros(iter_restart,numpy.float64)
944     else:     s=numpy.zeros(iter_restart,numpy.float64)
945        r=Msolve(b-Aprod(x))     g=numpy.zeros(iter_restart+1,numpy.float64)
       r_dot_r = bilinearform(r, r)  
       if r_dot_r<0: raise NegativeNorm,"negative norm."  
     
    h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)  
    c=numarray.zeros(iter_restart,numarray.Float64)  
    s=numarray.zeros(iter_restart,numarray.Float64)  
    g=numarray.zeros(iter_restart,numarray.Float64)  
946     v=[]     v=[]
947    
948       r_dot_r = bilinearform(r, r)
949       if r_dot_r<0: raise NegativeNorm,"negative norm."
950     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
951      
952     v.append(r/rho)     v.append(r/rho)
953     g[0]=rho     g[0]=rho
954    
955     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
956       while not (rho<=atol or iter==iter_restart):
957    
958      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
959    
960      p=Msolve(Aprod(v[iter]))          if P_R!=None:
961                p=Aprod(P_R(v[iter]))
962            else:
963            p=Aprod(v[iter])
964      v.append(p)      v.append(p)
965    
966      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
967    
968    # Modified Gram-Schmidt
969        for j in range(iter+1):
970          h[j,iter]=bilinearform(v[j],v[iter+1])
971          v[iter+1]-=h[j,iter]*v[j]
972    
973  # Modified Gram-Schmidt      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
974      for j in range(iter+1):      v_norm2=h[iter+1,iter]
       h[j][iter]=bilinearform(v[j],v[iter+1])    
       v[iter+1]-=h[j][iter]*v[j]  
         
     h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))  
     v_norm2=h[iter+1][iter]  
975    
976  # Reorthogonalize if needed  # Reorthogonalize if needed
977      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)
978       for j in range(iter+1):         for j in range(iter+1):
979          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
980              h[j][iter]=h[j][iter]+hr              h[j,iter]=h[j,iter]+hr
981              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
982    
983       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
984       h[iter+1][iter]=v_norm2       h[iter+1,iter]=v_norm2
985    
986  #   watch out for happy breakdown  #   watch out for happy breakdown
987          if not v_norm2 == 0:          if not v_norm2 == 0:
988           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
989    
990  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
991      if iter > 0 :      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
992          hhat=numarray.zeros(iter+1,numarray.Float64)      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
         for i in range(iter+1) : hhat[i]=h[i][iter]  
         hhat=__givapp(c[0:iter],s[0:iter],hhat);  
             for i in range(iter+1) : h[i][iter]=hhat[i]  
   
     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])  
993    
994      if mu!=0 :      if mu!=0 :
995          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
996          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
997          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
998          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
999          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1000                    g[iter]=gg[0]
1001                    g[iter+1]=gg[1]
1002  # Update the residual norm  # Update the residual norm
1003                  
1004          rho=abs(g[iter+1])          rho=abs(g[iter+1])
1005            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
1006      iter+=1      iter+=1
1007    
1008  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
1009  # It's time to compute x and leave.          # It's time to compute x and leave.
1010    
1011     if iter > 0 :     if verbose: print "GMRES: iteration stopped after %s step."%iter
1012       y=numarray.zeros(iter,numarray.Float64)         if iter > 0 :
1013       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y=numpy.zeros(iter,numpy.float64)
1014       if iter > 1 :         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
1015          i=iter-2         if iter > 1 :
1016            i=iter-2
1017          while i>=0 :          while i>=0 :
1018            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
1019            i=i-1            i=i-1
1020       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
1021       for i in range(iter-1):       for i in range(iter-1):
1022      xhat += v[i]*y[i]      xhat += v[i]*y[i]
1023     else : xhat=v[0]     else:
1024         xhat=v[0] * 0
1025     x += xhat     if P_R!=None:
1026     if iter<iter_restart-1:        x += P_R(xhat)
1027        stopped=True     else:
1028     else:        x += xhat
1029       if iter<iter_restart-1:
1030          stopped=True
1031       else:
1032        stopped=False        stopped=False
1033    
1034     return x,stopped     return x,stopped
1035    
1036  #################################################  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1037  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):      """
1038  #################################################      Solver for
1039      #  
1040      #  minres solves the system of linear equations Ax = b      *Ax=b*
1041      #  where A is a symmetric matrix (possibly indefinite or singular)  
1042      #  and b is a given vector.      with a symmetric and positive definite operator A (more details required!).
1043      #        It uses the minimum residual method (MINRES) with preconditioner M
1044      #  "A" may be a dense or sparse matrix (preferably sparse!)      providing an approximation of A.
1045      #  or the name of a function such that  
1046      #               y = A(x)      The iteration is terminated if
1047      #  returns the product y = Ax for any given vector x.  
1048      #      *|r| <= atol+rtol*|r0|*
1049      #  "M" defines a positive-definite preconditioner M = C C'.  
1050      #  "M" may be a dense or sparse matrix (preferably sparse!)      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1051      #  or the name of a function such that  
1052      #  solves the system My = x for any given vector x.      *|r| = sqrt( bilinearform(Msolve(r),r))*
1053      #  
1054      #      For details on the preconditioned conjugate gradient method see the book:
1055        
1056        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1057        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1058        C. Romine, and H. van der Vorst}.
1059    
1060        :param r: initial residual *r=b-Ax*. ``r`` is altered.
1061        :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1062        :param x: an initial guess for the solution
1063        :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1064        :param Aprod: returns the value Ax
1065        :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1066                     argument ``x``. The returned object needs to be of the same
1067                     type like argument ``r``.
1068        :param Msolve: solves Mx=r
1069        :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1070                      argument ``r``. The returned object needs to be of the same
1071                      type like argument ``x``.
1072        :param bilinearform: inner product ``<x,r>``
1073        :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1074                            type like argument ``x`` and ``r`` is. The returned value
1075                            is a ``float``.
1076        :param atol: absolute tolerance
1077        :type atol: non-negative ``float``
1078        :param rtol: relative tolerance
1079        :type rtol: non-negative ``float``
1080        :param iter_max: maximum number of iteration steps
1081        :type iter_max: ``int``
1082        :return: the solution approximation and the corresponding residual
1083        :rtype: ``tuple``
1084        :warning: ``r`` and ``x`` are altered.
1085        """
1086      #------------------------------------------------------------------      #------------------------------------------------------------------
1087      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
1088      # y  =  beta1 P' v1,  where  P = C**(-1).      # y  =  beta1 P' v1,  where  P = C**(-1).
1089      # v is really P' v1.      # v is really P' v1.
1090      #------------------------------------------------------------------      #------------------------------------------------------------------
1091      if x==None:      r1    = r
1092        x=0*b      y = Msolve(r)
1093      else:      beta1 = bilinearform(y,r)
       b += (-1)*Aprod(x)  
1094    
     r1    = b  
     y = Msolve(b)  
     beta1 = bilinearform(y,b)  
   
1095      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
1096    
1097      #  If b = 0 exactly, stop with x = 0.      #  If r = 0 exactly, stop with x
1098      if beta1==0: return x*0.      if beta1==0: return x
1099    
1100      if beta1> 0:      if beta1> 0: beta1  = math.sqrt(beta1)
       beta1  = math.sqrt(beta1)        
1101    
1102      #------------------------------------------------------------------      #------------------------------------------------------------------
1103      # Initialize quantities.      # Initialize quantities.
# Line 1008  def MINRES(b, Aprod, Msolve, bilinearfor Line 1117  def MINRES(b, Aprod, Msolve, bilinearfor
1117      ynorm2 = 0      ynorm2 = 0
1118      cs     = -1      cs     = -1
1119      sn     = 0      sn     = 0
1120      w      = b*0.      w      = r*0.
1121      w2     = b*0.      w2     = r*0.
1122      r2     = r1      r2     = r1
1123      eps    = 0.0001      eps    = 0.0001
1124    
1125      #---------------------------------------------------------------------      #---------------------------------------------------------------------
1126      # Main iteration loop.      # Main iteration loop.
1127      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1128      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1129    
1130      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
1131          iter    = iter  +  1          iter    = iter  +  1
# Line 1035  def MINRES(b, Aprod, Msolve, bilinearfor Line 1144  def MINRES(b, Aprod, Msolve, bilinearfor
1144          #-----------------------------------------------------------------          #-----------------------------------------------------------------
1145          s = 1/beta                 # Normalize previous vector (in y).          s = 1/beta                 # Normalize previous vector (in y).
1146          v = s*y                    # v = vk if P = I          v = s*y                    # v = vk if P = I
1147        
1148          y      = Aprod(v)          y      = Aprod(v)
1149        
1150          if iter >= 2:          if iter >= 2:
1151            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1152    
1153          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1154          y      += (- alfa/beta)*r2          y      += (- alfa/beta)*r2
1155          r1     = r2          r1     = r2
1156          r2     = y          r2     = y
1157          y = Msolve(r2)          y = Msolve(r2)
# Line 1052  def MINRES(b, Aprod, Msolve, bilinearfor Line 1161  def MINRES(b, Aprod, Msolve, bilinearfor
1161    
1162          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1163          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1164            
1165          if iter==1:                 # Initialize a few things.          if iter==1:                 # Initialize a few things.
1166            gmax   = abs( alfa )      # alpha1            gmax   = abs( alfa )      # alpha1
1167            gmin   = gmax             # alpha1            gmin   = gmax             # alpha1
# Line 1060  def MINRES(b, Aprod, Msolve, bilinearfor Line 1169  def MINRES(b, Aprod, Msolve, bilinearfor
1169          # Apply previous rotation Qk-1 to get          # Apply previous rotation Qk-1 to get
1170          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1171          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1172        
1173          oldeps = epsln          oldeps = epsln
1174          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1175          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
# Line 1070  def MINRES(b, Aprod, Msolve, bilinearfor Line 1179  def MINRES(b, Aprod, Msolve, bilinearfor
1179          # Compute the next plane rotation Qk          # Compute the next plane rotation Qk
1180    
1181          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1182          gamma  = max(gamma,eps)          gamma  = max(gamma,eps)
1183          cs     = gbar / gamma             # ck          cs     = gbar / gamma             # ck
1184          sn     = beta / gamma             # sk          sn     = beta / gamma             # sk
1185          phi    = cs * phibar              # phik          phi    = cs * phibar              # phik
# Line 1078  def MINRES(b, Aprod, Msolve, bilinearfor Line 1187  def MINRES(b, Aprod, Msolve, bilinearfor
1187    
1188          # Update  x.          # Update  x.
1189    
1190          denom = 1/gamma          denom = 1/gamma
1191          w1    = w2          w1    = w2
1192          w2    = w          w2    = w
1193          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1194          x     +=  phi*w          x     +=  phi*w
1195    
# Line 1095  def MINRES(b, Aprod, Msolve, bilinearfor Line 1204  def MINRES(b, Aprod, Msolve, bilinearfor
1204    
1205          # Estimate various norms and test for convergence.          # Estimate various norms and test for convergence.
1206    
1207          Anorm  = math.sqrt( tnorm2 )          Anorm  = math.sqrt( tnorm2 )
1208          ynorm  = math.sqrt( ynorm2 )          ynorm  = math.sqrt( ynorm2 )
1209    
1210          rnorm  = phibar          rnorm  = phibar
1211    
1212      return x      return x
1213    
1214  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1215      """
1216      Solver for
1217    
1218  # TFQMR solver for linear systems    *Ax=b*
1219  #  
1220  #    with a general operator A (more details required!).
1221  # initialization    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
 #  
   errtol = math.sqrt(bilinearform(b,b))  
   norm_b=errtol  
   kmax  = iter_max  
   error = []  
   
   if math.sqrt(bilinearform(x,x)) != 0.0:  
     r = b - Aprod(x)  
   else:  
     r = b  
1222    
1223    r=Msolve(r)    The iteration is terminated if
1224    
1225      *|r| <= atol+rtol*|r0|*
1226    
1227      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1228    
1229      *|r| = sqrt( bilinearform(r,r))*
1230    
1231      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1232      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1233      :param x: an initial guess for the solution
1234      :type x: same like ``r``
1235      :param Aprod: returns the value Ax
1236      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1237                   argument ``x``. The returned object needs to be of the same type
1238                   like argument ``r``.
1239      :param bilinearform: inner product ``<x,r>``
1240      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1241                          type like argument ``x`` and ``r``. The returned value is
1242                          a ``float``.
1243      :param atol: absolute tolerance
1244      :type atol: non-negative ``float``
1245      :param rtol: relative tolerance
1246      :type rtol: non-negative ``float``
1247      :param iter_max: maximum number of iteration steps
1248      :type iter_max: ``int``
1249      :rtype: ``tuple``
1250      :warning: ``r`` and ``x`` are altered.
1251      """
1252    u1=0    u1=0
1253    u2=0    u2=0
1254    y1=0    y1=0
1255    y2=0    y2=0
1256    
1257    w = r    w = r
1258    y1 = r    y1 = r
1259    iter = 0    iter = 0
1260    d = 0    d = 0
1261        v = Aprod(y1)
   v = Msolve(Aprod(y1))  
1262    u1 = v    u1 = v
1263      
1264    theta = 0.0;    theta = 0.0;
1265    eta = 0.0;    eta = 0.0;
1266    tau = math.sqrt(bilinearform(r,r))    rho=bilinearform(r,r)
1267    error = [ error, tau ]    if rho < 0: raise NegativeNorm,"negative norm."
1268    rho = tau * tau    tau = math.sqrt(rho)
1269    m=1    norm_r0=tau
1270  #    while tau>atol+rtol*norm_r0:
 #  TFQMR iteration  
 #  
 #  while ( iter < kmax-1 ):  
     
   while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):  
1271      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
1272    
1273      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1274        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
     if ( sigma == 0.0 ):  
       raise 'TFQMR breakdown, sigma=0'  
       
1275    
1276      alpha = rho / sigma      alpha = rho / sigma
1277    
# Line 1162  def TFQMR(b, Aprod, Msolve, bilinearform Line 1281  def TFQMR(b, Aprod, Msolve, bilinearform
1281  #  #
1282        if ( j == 1 ):        if ( j == 1 ):
1283          y2 = y1 - alpha * v          y2 = y1 - alpha * v
1284          u2 = Msolve(Aprod(y2))          u2 = Aprod(y2)
1285    
1286        m = 2 * (iter+1) - 2 + (j+1)        m = 2 * (iter+1) - 2 + (j+1)
1287        if j==0:        if j==0:
1288           w = w - alpha * u1           w = w - alpha * u1
1289           d = y1 + ( theta * theta * eta / alpha ) * d           d = y1 + ( theta * theta * eta / alpha ) * d
1290        if j==1:        if j==1:
# Line 1180  def TFQMR(b, Aprod, Msolve, bilinearform Line 1299  def TFQMR(b, Aprod, Msolve, bilinearform
1299  #  #
1300  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1301  #  #
1302       # if ( tau * math.sqrt ( m + 1 ) <= errtol ):      if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
      #   error = [ error, tau ]  
      #   total_iters = iter  
      #   break  
         
   
     if ( rho == 0.0 ):  
       raise 'TFQMR breakdown, rho=0'  
       
1303    
1304      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1305      beta = rhon / rho;      beta = rhon / rho;
1306      rho = rhon;      rho = rhon;
1307      y1 = w + beta * y2;      y1 = w + beta * y2;
1308      u1 = Msolve(Aprod(y1))      u1 = Aprod(y1)
1309      v = u1 + beta * ( u2 + beta * v )      v = u1 + beta * ( u2 + beta * v )
1310      error = [ error, tau ]  
1311      total_iters = iter      iter += 1
       
     iter = iter + 1  
1312    
1313    return x    return x
1314    
# Line 1208  def TFQMR(b, Aprod, Msolve, bilinearform Line 1317  def TFQMR(b, Aprod, Msolve, bilinearform
1317    
1318  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1319     """     """
1320     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.     Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1321       ArithmeticTuple and ``a`` is a float.
1322    
1323     example of usage:     Example of usage::
1324    
1325     from esys.escript import Data         from esys.escript import Data
1326     from numarray import array         from numpy import array
1327     a=Data(...)         a=Data(...)
1328     b=array([1.,4.])         b=array([1.,4.])
1329     x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1330     y=5.*x         y=5.*x
1331    
1332     """     """
1333     def __init__(self,*args):     def __init__(self,*args):
1334         """         """
1335         initialize object with elements args.         Initializes object with elements ``args``.
1336    
1337         @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)         :param args: tuple of objects that support inplace add (x+=y) and
1338                        scaling (x=a*y)
1339         """         """
1340         self.__items=list(args)         self.__items=list(args)
1341    
1342     def __len__(self):     def __len__(self):
1343         """         """
1344         number of items         Returns the number of items.
1345    
1346         @return: number of items         :return: number of items
1347         @rtype: C{int}         :rtype: ``int``
1348         """         """
1349         return len(self.__items)         return len(self.__items)
1350    
1351     def __getitem__(self,index):     def __getitem__(self,index):
1352         """         """
1353         get an item         Returns item at specified position.
1354    
1355         @param index: item to be returned         :param index: index of item to be returned
1356         @type index: C{int}         :type index: ``int``
1357         @return: item with index C{index}         :return: item with index ``index``
1358         """         """
1359         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1360    
1361     def __mul__(self,other):     def __mul__(self,other):
1362         """         """
1363         scaling from the right         Scales by ``other`` from the right.
1364    
1365         @param other: scaling factor         :param other: scaling factor
1366         @type other: C{float}         :type other: ``float``
1367         @return: itemwise self*other         :return: itemwise self*other
1368         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1369         """         """
1370         out=[]         out=[]
1371         try:           try:
1372             l=len(other)             l=len(other)
1373             if l!=len(self):             if l!=len(self):
1374                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1375             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1376         except TypeError:         except TypeError:
1377         for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
1378         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1379    
1380     def __rmul__(self,other):     def __rmul__(self,other):
1381         """         """
1382         scaling from the left         Scales by ``other`` from the left.
1383    
1384         @param other: scaling factor         :param other: scaling factor
1385         @type other: C{float}         :type other: ``float``
1386         @return: itemwise other*self         :return: itemwise other*self
1387         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1388         """         """
1389         out=[]         out=[]
1390         try:           try:
1391             l=len(other)             l=len(other)
1392             if l!=len(self):             if l!=len(self):
1393                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1394             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1395         except TypeError:         except TypeError:
1396         for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
1397         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1398    
1399     def __div__(self,other):     def __div__(self,other):
1400         """         """
1401         dividing from the right         Scales by (1/``other``) from the right.
1402    
1403         @param other: scaling factor         :param other: scaling factor
1404         @type other: C{float}         :type other: ``float``
1405         @return: itemwise self/other         :return: itemwise self/other
1406         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1407         """         """
1408         return self*(1/other)         return self*(1/other)
1409    
1410     def __rdiv__(self,other):     def __rdiv__(self,other):
1411         """         """
1412         dividing from the left         Scales by (1/``other``) from the left.
1413    
1414         @param other: scaling factor         :param other: scaling factor
1415         @type other: C{float}         :type other: ``float``
1416         @return: itemwise other/self         :return: itemwise other/self
1417         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1418         """         """
1419         out=[]         out=[]
1420         try:           try:
1421             l=len(other)             l=len(other)
1422             if l!=len(self):             if l!=len(self):
1423                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1424             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1425         except TypeError:         except TypeError:
1426         for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
1427         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1428      
1429     def __iadd__(self,other):     def __iadd__(self,other):
1430         """         """
1431         in-place add of other to self         Inplace addition of ``other`` to self.
1432    
1433         @param other: increment         :param other: increment
1434         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1435         """         """
1436         if len(self) != len(other):         if len(self) != len(other):
1437             raise ValueError,"tuple length must match."             raise ValueError,"tuple lengths must match."
1438         for i in range(len(self)):         for i in range(len(self)):
1439             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1440         return self         return self
1441    
1442     def __add__(self,other):     def __add__(self,other):
1443         """         """
1444         add other to self         Adds ``other`` to self.
1445    
1446         @param other: increment         :param other: increment
1447         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1448         """         """
1449         out=[]         out=[]
1450         try:           try:
1451             l=len(other)             l=len(other)
1452             if l!=len(self):             if l!=len(self):
1453                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1454             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1455         except TypeError:         except TypeError:
1456         for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
1457         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1458    
1459     def __sub__(self,other):     def __sub__(self,other):
1460         """         """
1461         subtract other from self         Subtracts ``other`` from self.
1462    
1463         @param other: increment         :param other: decrement
1464         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1465         """         """
1466         out=[]         out=[]
1467         try:           try:
1468             l=len(other)             l=len(other)
1469             if l!=len(self):             if l!=len(self):
1470                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1471             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1472         except TypeError:         except TypeError:
1473         for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
1474         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1475      
1476     def __isub__(self,other):     def __isub__(self,other):
1477         """         """
1478         subtract other from self         Inplace subtraction of ``other`` from self.
1479    
1480         @param other: increment         :param other: decrement
1481         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1482         """         """
1483         if len(self) != len(other):         if len(self) != len(other):
1484             raise ValueError,"tuple length must match."             raise ValueError,"tuple length must match."
# Line 1377  class ArithmeticTuple(object): Line 1488  class ArithmeticTuple(object):
1488    
1489     def __neg__(self):     def __neg__(self):
1490         """         """
1491         negate         Negates values.
   
1492         """         """
1493         out=[]         out=[]
1494         for i in range(len(self)):         for i in range(len(self)):
# Line 1388  class ArithmeticTuple(object): Line 1498  class ArithmeticTuple(object):
1498    
1499  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1500        """        """
1501        This provides a framwork for solving linear homogeneous saddle point problem of the form        This class provides a framework for solving linear homogeneous saddle
1502          point problems of the form::
              Av+B^*p=f  
              Bv    =0  
1503    
1504        for the unknowns v and p and given operators A and B and given right hand side f.            *Av+B^*p=f*
1505        B^* is the adjoint operator of B is the given inner product.            *Bv     =0*
1506    
1507          for the unknowns *v* and *p* and given operators *A* and *B* and
1508          given right hand side *f*. *B^** is the adjoint operator of *B*.
1509          *A* may depend weakly on *v* and *p*.
1510        """        """
1511        def __init__(self,**kwargs):        def __init__(self, **kwargs):
1512        """
1513        initializes the saddle point problem
1514        """
1515            self.resetControlParameters()
1516          self.setTolerance()          self.setTolerance()
1517          self.setToleranceReductionFactor()          self.setAbsoluteTolerance()
1518          def resetControlParameters(self, K_p=1., K_v=1., rtol_max=0.01, rtol_min = 1.e-7, chi_max=0.5, reduction_factor=0.3, theta = 0.1):
       def initialize(self):  
         """  
         initialize the problem (overwrite)  
         """  
         pass  
       def B(self,v):  
1519           """           """
1520           returns Bv (overwrite)           sets a control parameter
          @rtype: equal to the type of p  
1521    
1522           @note: boundary conditions on p should be zero!           :param K_p: initial value for constant to adjust pressure tolerance
1523             :type K_p: ``float``
1524             :param K_v: initial value for constant to adjust velocity tolerance
1525             :type K_v: ``float``
1526             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1527             :type rtol_max: ``float``
1528             :param chi_max: maximum tolerable converegence rate.
1529             :type chi_max: ``float``
1530             :param reduction_factor: reduction factor for adjustment factors.
1531             :type reduction_factor: ``float``
1532           """           """
1533           pass           self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1534    
1535        def inner(self,p0,p1):        def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1536           """           """
1537           returns inner product of two element p0 and p1  (overwrite)           sets a control parameter
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
1538    
          @rtype: equal to the type of p  
          """  
          pass  
1539    
1540        def solve_A(self,u,p):           :param K_p: initial value for constant to adjust pressure tolerance
1541             :type K_p: ``float``
1542             :param K_v: initial value for constant to adjust velocity tolerance
1543             :type K_v: ``float``
1544             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1545             :type rtol_max: ``float``
1546             :param chi_max: maximum tolerable converegence rate.
1547             :type chi_max: ``float``
1548             :type reduction_factor: ``float``
1549           """           """
1550           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           if not K_p == None:
1551                if K_p<1:
1552                   raise ValueError,"K_p need to be greater or equal to 1."
1553             else:
1554                K_p=self.__K_p
1555    
1556             if not K_v == None:
1557                if K_v<1:
1558                   raise ValueError,"K_v need to be greater or equal to 1."
1559             else:
1560                K_v=self.__K_v
1561    
1562             if not rtol_max == None:
1563                if rtol_max<=0 or rtol_max>=1:
1564                   raise ValueError,"rtol_max needs to be positive and less than 1."
1565             else:
1566                rtol_max=self.__rtol_max
1567    
1568             if not rtol_min == None:
1569                if rtol_min<=0 or rtol_min>=1:
1570                   raise ValueError,"rtol_min needs to be positive and less than 1."
1571             else:
1572                rtol_min=self.__rtol_min
1573    
1574             if not chi_max == None:
1575                if chi_max<=0 or chi_max>=1:
1576                   raise ValueError,"chi_max needs to be positive and less than 1."
1577             else:
1578                chi_max = self.__chi_max
1579    
1580             if not reduction_factor == None:
1581                if reduction_factor<=0 or reduction_factor>1:
1582                   raise ValueError,"reduction_factor need to be between zero and one."
1583             else:
1584                reduction_factor=self.__reduction_factor
1585    
1586             if not theta == None:
1587                if theta<=0 or theta>1:
1588                   raise ValueError,"theta need to be between zero and one."
1589             else:
1590                theta=self.__theta
1591    
1592             if rtol_min>=rtol_max:
1593                 raise ValueError,"rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min)
1594             self.__chi_max = chi_max
1595             self.__rtol_max = rtol_max
1596             self.__K_p = K_p
1597             self.__K_v = K_v
1598             self.__reduction_factor = reduction_factor
1599             self.__theta = theta
1600             self.__rtol_min=rtol_min
1601    
1602           @rtype: equal to the type of v        #=============================================================
1603           @note: boundary conditions on v should be zero!        def inner_pBv(self,p,Bv):
1604           """           """
1605           pass           Returns inner product of element p and Bv (overwrite).
1606    
1607        def solve_prec(self,p):           :param p: a pressure increment
1608             :param Bv: a residual
1609             :return: inner product of element p and Bv
1610             :rtype: ``float``
1611             :note: used if PCG is applied.
1612           """           """
1613           provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)           raise NotImplementedError,"no inner product for p and Bv implemented."
1614    
1615           @rtype: equal to the type of p        def inner_p(self,p0,p1):
1616           """           """
1617           pass           Returns inner product of p0 and p1 (overwrite).
1618    
1619        def stoppingcriterium(self,Bv,v,p):           :param p0: a pressure
1620             :param p1: a pressure
1621             :return: inner product of p0 and p1
1622             :rtype: ``float``
1623           """           """
1624           returns a True if iteration is terminated. (overwrite)           raise NotImplementedError,"no inner product for p implemented."
1625      
1626          def norm_v(self,v):
1627             """
1628             Returns the norm of v (overwrite).
1629    
1630           @rtype: C{bool}           :param v: a velovity
1631             :return: norm of v
1632             :rtype: non-negative ``float``
1633             """
1634             raise NotImplementedError,"no norm of v implemented."
1635          def getDV(self, p, v, tol):
1636           """           """
1637           pass           return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
               
       def __inner(self,p,r):  
          return self.inner(p,r[1])  
1638    
1639        def __inner_p(self,p1,p2):           :param p: pressure
1640           return self.inner(p1,p2)           :param v: pressure
1641                   :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1642        def __inner_a(self,a1,a2):           :note: Only *A* may depend on *v* and *p*
1643           return self.inner_a(a1,a2)           """
1644             raise NotImplementedError,"no dv calculation implemented."
1645    
1646        def __inner_a1(self,a1,a2):          
1647           return self.inner(a1[1],a2[1])        def Bv(self,v, tol):
1648            """
1649            Returns Bv with accuracy `tol` (overwrite)
1650    
1651        def __stoppingcriterium(self,norm_r,r,p):          :rtype: equal to the type of p
1652            return self.stoppingcriterium(r[1],r[0],p)          :note: boundary conditions on p should be zero!
1653            """
1654            raise NotImplementedError, "no operator B implemented."
1655    
1656        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):        def norm_Bv(self,Bv):
1657            return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)          """
1658            Returns the norm of Bv (overwrite).
1659    
1660        def setTolerance(self,tolerance=1.e-8):          :rtype: equal to the type of p
1661                self.__tol=tolerance          :note: boundary conditions on p should be zero!
1662        def getTolerance(self):          """
1663                return self.__tol          raise NotImplementedError, "no norm of Bv implemented."
       def setToleranceReductionFactor(self,reduction=0.01):  
               self.__reduction=reduction  
       def getSubProblemTolerance(self):  
               return self.__reduction*self.getTolerance()  
   
       def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):  
               """  
               solves the saddle point problem using initial guesses v and p.  
   
               @param max_iter: maximum number of iteration steps.  
               """  
               self.verbose=verbose  
               self.show_details=show_details and self.verbose  
   
               # assume p is known: then v=A^-1(f-B^*p)  
               # which leads to BA^-1B^*p = BA^-1f    
   
           # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)        
           self.__z=v+self.solve_A(v,p*0)  
               Bz=self.B(self.__z)  
               #  
           #   solve BA^-1B^*p = Bz  
               #  
               #  
               #  
               self.iter=0  
           if solver=='GMRES':        
                 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter  
                 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='TFQMR':        
                 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter  
                 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='MINRES':        
                 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter  
                 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
                 
           if solver=='GMRESC':        
                 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter  
                 p0=self.solve_prec1(Bz)  
             #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))  
                 #p-=alfa  
                 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=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())  
   
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
             p=x[1]  
         u=v+self.solve_A(v,p)        
         #p=x[1]  
         #u=x[0]  
   
               if solver=='PCG':  
                 #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
                 #  
                 #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)  
                 #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)  
                 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter  
                 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)  
             u=r[0]    
                 # print "DIFF=",util.integrate(p)  
   
               # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)  
   
           return u,p  
   
       def __Msolve(self,r):  
           return self.solve_prec(r[1])  
   
       def __Msolve2(self,r):  
           return self.solve_prec(r*1)  
   
       def __Mempty(self,r):  
           return r  
   
   
       def __Aprod(self,p):  
           # return BA^-1B*p  
           #solve Av =B^*p as Av =f-Az-B^*(-p)  
           v=self.solve_A(self.__z,-p)  
           return ArithmeticTuple(v, self.B(v))  
   
       def __Anext(self,x):  
           # return next v,p  
           #solve Av =-B^*p as Av =f-Az-B^*p  
   
       pc=x[1]  
           v=self.solve_A(self.__z,-pc)  
       p=self.solve_prec1(self.B(v))  
   
           return ArithmeticTuple(v,p)  
   
   
       def __Aprod2(self,p):  
           # return BA^-1B*p  
           #solve Av =B^*p as Av =f-Az-B^*(-p)  
       v=self.solve_A(self.__z,-p)  
           return self.B(v)  
   
       def __Aprod_Newton(self,p):  
           # return BA^-1B*p - Bz  
           #solve Av =-B^*p as Av =f-Az-B^*p  
       v=self.solve_A(self.__z,-p)  
           return self.B(v-self.__z)  
   
       def __Aprod_Newton2(self,x):  
           # return BA^-1B*p - Bz  
           #solve Av =-B^*p as Av =f-Az-B^*p  
           pc=x[1]  
       v=self.solve_A(self.__z,-pc)  
           p=self.solve_prec1(self.B(v-self.__z))  
           return ArithmeticTuple(v,p)  
1664    
1665          def solve_AinvBt(self,dp, tol):
1666             """
1667             Solves *A dv=B^*dp* with accuracy `tol`
1668    
1669  def MaskFromBoundaryTag(domain,*tags):           :param dp: a pressure increment
1670     """           :return: the solution of *A dv=B^*dp*
1671     creates a mask on the Solution(domain) function space which one for samples           :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1672     that touch the boundary tagged by tags.           """
1673             raise NotImplementedError,"no operator A implemented."
1674    
1675     usage: m=MaskFromBoundaryTag(domain,"left", "right")        def solve_prec(self,Bv, tol):
1676             """
1677             Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
1678    
1679     @param domain: a given domain           :rtype: equal to the type of p
1680     @type domain: L{escript.Domain}           :note: boundary conditions on p should be zero!
1681     @param tags: boundray tags           """
1682     @type tags: C{str}           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1683     @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.        #=============================================================
1684     @rtype: L{escript.Data} of rank 0        def __Aprod_PCG(self,dp):
1685     """            dv=self.solve_AinvBt(dp, self.__subtol)
1686     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1687     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))  
1688     for t in tags: d.setTaggedValue(t,1.)        def __inner_PCG(self,p,r):
1689     pde.setValue(y=d)           return self.inner_pBv(p,r[1])
1690     return util.whereNonZero(pde.getRightHandSide())  
1691  #============================================================================================================================================        def __Msolve_PCG(self,r):
1692  class SaddlePointProblem(object):            return self.solve_prec(r[1], self.__subtol)
1693     """        #=============================================================
1694     This implements a solver for a saddlepoint problem        def __Aprod_GMRES(self,p):
1695              return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1696          def __inner_GMRES(self,p0,p1):
1697             return self.inner_p(p0,p1)
1698    
1699          #=============================================================
1700          def norm_p(self,p):
1701              """
1702              calculates the norm of ``p``
1703              
1704              :param p: a pressure
1705              :return: the norm of ``p`` using the inner product for pressure
1706              :rtype: ``float``
1707              """
1708              f=self.inner_p(p,p)
1709              if f<0: raise ValueError,"negative pressure norm."
1710              return math.sqrt(f)
1711          
1712          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1713             """
1714             Solves the saddle point problem using initial guesses v and p.
1715    
1716     M{f(u,p)=0}           :param v: initial guess for velocity
1717     M{g(u)=0}           :param p: initial guess for pressure
1718             :type v: `Data`
1719             :type p: `Data`
1720             :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1721             :param max_iter: maximum number of iteration steps per correction
1722                              attempt
1723             :param verbose: if True, shows information on the progress of the
1724                             saddlepoint problem solver.
1725             :param iter_restart: restart the iteration after ``iter_restart`` steps
1726                                  (only used if useUzaw=False)
1727             :type usePCG: ``bool``
1728             :type max_iter: ``int``
1729             :type verbose: ``bool``
1730             :type iter_restart: ``int``
1731             :rtype: ``tuple`` of `Data` objects
1732             :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1733             """
1734             return self._solve(v=v,p=p,max_iter=max_iter,verbose=verbose, usePCG=usePCG, iter_restart=iter_restart, max_correction_steps=max_correction_steps)
1735    
1736     for u and p. The problem is solved with an inexact Uszawa scheme for p:        def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1737             """
1738             see `_solve` method.
1739             """
1740             self.verbose=verbose
1741             rtol=self.getTolerance()
1742             atol=self.getAbsoluteTolerance()
1743    
1744             K_p=self.__K_p
1745             K_v=self.__K_v
1746             correction_step=0
1747             converged=False
1748             chi=None
1749             eps=None
1750    
1751             if self.verbose: print "HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)
1752             while not converged:
1753    
1754                 # get tolerance for velecity increment:
1755                 if chi == None:
1756                    rtol_v=self.__rtol_max
1757                 else:
1758                    rtol_v=min(chi/K_v,self.__rtol_max)
1759                 rtol_v=max(rtol_v, self.__rtol_min)
1760                 if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)
1761                 # get velocity increment:
1762                 dv1=self.getDV(p,v,rtol_v)
1763                 v1=v+dv1
1764                 Bv1=self.Bv(v1, rtol_v)
1765                 norm_Bv1=self.norm_Bv(Bv1)
1766                 norm_dv1=self.norm_v(dv1)
1767                 if self.verbose: print "HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)
1768                 if norm_dv1*self.__theta < norm_Bv1:
1769                    # get tolerance for pressure increment:
1770                    large_Bv1=True
1771                    if chi == None or eps == None:
1772                       rtol_p=self.__rtol_max
1773                    else:
1774                       rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1775                    self.__subtol=max(rtol_p**2, self.__rtol_min)
1776                    if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)
1777                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1778                    if usePCG:
1779                        dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)
1780                        v2=r[0]
1781                        Bv2=r[1]
1782                    else:
1783                        # don't use!!!!
1784                        dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1785                        dv2=self.solve_AinvBt(dp, self.__subtol)
1786                        v2=v1-dv2
1787                        Bv2=self.Bv(v2, self.__subtol)
1788                    p2=p+dp
1789                 else:
1790                    large_Bv1=False
1791                    v2=v1
1792                    p2=p
1793                 # update business:
1794                 norm_dv2=self.norm_v(v2-v)
1795                 norm_v2=self.norm_v(v2)
1796                 if self.verbose: print "HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))
1797                 eps, eps_old = max(norm_Bv1, norm_dv2), eps
1798                 if eps_old == None:
1799                      chi, chi_old = None, chi
1800                 else:
1801                      chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1802                 if eps != None:
1803                     if chi !=None:
1804                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)
1805                     else:
1806                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)
1807                 if eps <= rtol*norm_v2+atol :
1808                     converged = True
1809                 else:
1810                     if correction_step>=max_correction_steps:
1811                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1812                     if chi_old!=None:
1813                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1814                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1815                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)
1816                 correction_step+=1
1817                 v,p =v2, p2
1818             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1819         return v,p
1820          #========================================================================
1821          def setTolerance(self,tolerance=1.e-4):
1822             """
1823             Sets the relative tolerance for (v,p).
1824    
1825     M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}           :param tolerance: tolerance to be used
1826     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}           :type tolerance: non-negative ``float``
1827             """
1828             if tolerance<0:
1829                 raise ValueError,"tolerance must be positive."
1830             self.__rtol=tolerance
1831    
1832     where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of        def getTolerance(self):
1833     A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'           """
1834     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays           Returns the relative tolerance.
    in fact the role of a preconditioner.  
    """  
    def __init__(self,verbose=False,*args):  
        """  
        initializes the problem  
1835    
1836         @param verbose: switches on the printing out some information           :return: relative tolerance
1837         @type verbose: C{bool}           :rtype: ``float``
1838         @note: this method may be overwritten by a particular saddle point problem           """
1839         """           return self.__rtol
        print "SaddlePointProblem should not be used anymore!"  
        if not isinstance(verbose,bool):  
             raise TypeError("verbose needs to be of type bool.")  
        self.__verbose=verbose  
        self.relaxation=1.  
        DeprecationWarning("SaddlePointProblem should not be used anymore.")  
1840    
1841     def trace(self,text):        def setAbsoluteTolerance(self,tolerance=0.):
1842         """           """
1843         prints text if verbose has been set           Sets the absolute tolerance.
1844    
1845         @param text: a text message           :param tolerance: tolerance to be used
1846         @type text: C{str}           :type tolerance: non-negative ``float``
1847         """           """
1848         if self.__verbose: print "%s: %s"%(str(self),text)           if tolerance<0:
1849                 raise ValueError,"tolerance must be non-negative."
1850             self.__atol=tolerance
1851    
1852     def solve_f(self,u,p,tol=1.e-8):        def getAbsoluteTolerance(self):
1853         """           """
1854         solves           Returns the absolute tolerance.
1855    
1856         A_f du = f(u,p)           :return: absolute tolerance
1857             :rtype: ``float``
1858             """
1859             return self.__atol
1860    
        with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.  
1861    
1862         @param u: current approximation of u  def MaskFromBoundaryTag(domain,*tags):
1863         @type u: L{escript.Data}     """
1864         @param p: current approximation of p     Creates a mask on the Solution(domain) function space where the value is
1865         @type p: L{escript.Data}     one for samples that touch the boundary tagged by tags.
        @param tol: tolerance expected for du  
        @type tol: C{float}  
        @return: increment du  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point problem  
        """  
        pass  
1866    
1867     def solve_g(self,u,tol=1.e-8):     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
        """  
        solves  
1868    
1869         Q_g dp = g(u)     :param domain: domain to be used
1870       :type domain: `escript.Domain`
1871       :param tags: boundary tags
1872       :type tags: ``str``
1873       :return: a mask which marks samples that are touching the boundary tagged
1874                by any of the given tags
1875       :rtype: `escript.Data` of rank 0
1876       """
1877       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1878       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1879       for t in tags: d.setTaggedValue(t,1.)
1880       pde.setValue(y=d)
1881       return util.whereNonZero(pde.getRightHandSide())
1882    
1883         with Q_g is a preconditioner for A_g A_f^{-1} A_g with  A_g is the jacobiean of g with respect to p.  def MaskFromTag(domain,*tags):
1884       """
1885       Creates a mask on the Solution(domain) function space where the value is
1886       one for samples that touch regions tagged by tags.
1887    
1888         @param u: current approximation of u     Usage: m=MaskFromTag(domain, "ham")
        @type u: L{escript.Data}  
        @param tol: tolerance expected for dp  
        @type tol: C{float}  
        @return: increment dp  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point problem  
        """  
        pass  
1889    
1890     def inner(self,p0,p1):     :param domain: domain to be used
1891         """     :type domain: `escript.Domain`
1892         inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)     :param tags: boundary tags
1893         @return: inner product of p0 and p1     :type tags: ``str``
1894         @rtype: C{float}     :return: a mask which marks samples that are touching the boundary tagged
1895         """              by any of the given tags
1896         pass     :rtype: `escript.Data` of rank 0
1897       """
1898       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1899       d=escript.Scalar(0.,escript.Function(domain))
1900       for t in tags: d.setTaggedValue(t,1.)
1901       pde.setValue(Y=d)
1902       return util.whereNonZero(pde.getRightHandSide())
1903    
    subiter_max=3  
    def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):  
         """  
         runs the solver  
1904    
         @param u0: initial guess for C{u}  
         @type u0: L{esys.escript.Data}  
         @param p0: initial guess for C{p}  
         @type p0: L{esys.escript.Data}  
         @param tolerance: tolerance for relative error in C{u} and C{p}  
         @type tolerance: positive C{float}  
         @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}  
         @type tolerance_u: positive C{float}  
         @param iter_max: maximum number of iteration steps.  
         @type iter_max: C{int}  
         @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the  
                                    relaxation factor. If C{accepted_reduction=None} no backtracking is used.  
         @type accepted_reduction: positive C{float} or C{None}  
         @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.  
         @type relaxation: C{float} or C{None}  
         """  
         tol=1.e-2  
         if tolerance_u==None: tolerance_u=tolerance  
         if not relaxation==None: self.relaxation=relaxation  
         if accepted_reduction ==None:  
               angle_limit=0.  
         elif accepted_reduction>=1.:  
               angle_limit=0.  
         else:  
               angle_limit=util.sqrt(1-accepted_reduction**2)  
         self.iter=0  
         u=u0  
         p=p0  
         #  
         #   initialize things:  
         #  
         converged=False  
         #  
         #  start loop:  
         #  
         #  initial search direction is g  
         #  
         while not converged :  
             if self.iter>iter_max:  
                 raise ArithmeticError("no convergence after %s steps."%self.iter)  
             f_new=self.solve_f(u,p,tol)  
             norm_f_new = util.Lsup(f_new)  
             u_new=u-f_new  
             g_new=self.solve_g(u_new,tol)  
             self.iter+=1  
             norm_g_new = util.sqrt(self.inner(g_new,g_new))  
             if norm_f_new==0. and norm_g_new==0.: return u, p  
             if self.iter>1 and not accepted_reduction==None:  
                #  
                #   did we manage to reduce the norm of G? I  
                #   if not we start a backtracking procedure  
                #  
                # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g  
                if norm_g_new > accepted_reduction * norm_g:  
                   sub_iter=0  
                   s=self.relaxation  
                   d=g  
                   g_last=g  
                   self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))  
                   while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:  
                      dg= g_new-g_last  
                      norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)  
                      rad=self.inner(g_new,dg)/self.relaxation  
                      # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit  
                      # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g  
                      if abs(rad) < norm_dg*norm_g_new * angle_limit:  
                          if sub_iter>0: self.trace("    no further improvements expected from backtracking.")  
                          break  
                      r=self.relaxation  
                      self.relaxation= - rad/norm_dg**2  
                      s+=self.relaxation  
                      #####  
                      # a=g_new+self.relaxation*dg/r  
                      # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation  
                      #####  
                      g_last=g_new  
                      p+=self.relaxation*d  
                      f_new=self.solve_f(u,p,tol)  
                      u_new=u-f_new  
                      g_new=self.solve_g(u_new,tol)  
                      self.iter+=1  
                      norm_f_new = util.Lsup(f_new)  
                      norm_g_new = util.sqrt(self.inner(g_new,g_new))  
                      # print "   ",sub_iter," new g norm",norm_g_new  
                      self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))  
                      #  
                      #   can we expect reduction of g?  
                      #  
                      # u_last=u_new  
                      sub_iter+=1  
                   self.relaxation=s  
             #  
             #  check for convergence:  
             #  
             norm_u_new = util.Lsup(u_new)  
             p_new=p+self.relaxation*g_new  
             norm_p_new = util.sqrt(self.inner(p_new,p_new))  
             self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))  
   
             if self.iter>1:  
                dg2=g_new-g  
                df2=f_new-f  
                norm_dg2=util.sqrt(self.inner(dg2,dg2))  
                norm_df2=util.Lsup(df2)  
                # print norm_g_new, norm_g, norm_dg, norm_p, tolerance  
                tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new  
                tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new  
                if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:  
                    converged=True  
             f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new  
         self.trace("convergence after %s steps."%self.iter)  
         return u,p  

Legend:
Removed from v.1956  
changed lines
  Added in v.2881

  ViewVC Help
Powered by ViewVC 1.1.26