/[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 2719 by gross, Wed Oct 14 06:38:03 2009 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 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-2009 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 __call__(self, input_data):    def __call__(self, input_data):
157      """      """
158      Projects input_data onto a continuous function      Projects ``input_data`` onto a continuous function.
159    
160      @param input_data: The input_data to be projected.      :param input_data: the data to be projected
161      """      """
162      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
163      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 190  class Projector:
190    
191  class NoPDE:  class NoPDE:
192       """       """
193       solves the following problem for u:       Solves the following problem for u:
194    
195       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
196    
197       with constraint       with constraint
198    
199       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
200    
201       where D, Y, r and q are given functions of rank 1.       where *D*, *Y*, *r* and *q* are given functions of rank 1.
202    
203       In the case of scalars this takes the form       In the case of scalars this takes the form
204    
205       M{D*u=Y}       *D*u=Y*
206    
207       with constraint       with constraint
208    
209       M{u=r}  where M{q>0}       *u=r* where *q>0*
210    
211       where D, Y, r and q are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
212    
213       The constraint is overwriting any other condition.       The constraint overwrites any other condition.
214    
215       @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
216              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
217              thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.              in `Solution` or `ReducedSolution`.
218       """       """
219         # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
220         # this.
221       def __init__(self,domain,D=None,Y=None,q=None,r=None):       def __init__(self,domain,D=None,Y=None,q=None,r=None):
222           """           """
223           initialize the problem           Initializes the problem.
224    
225           @param domain: domain of the PDE.           :param domain: domain of the PDE
226           @type domain: L{Domain}           :type domain: `Domain`
227           @param D: coefficient of the solution.           :param D: coefficient of the solution
228           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
229           @param Y: right hand side           :param Y: right hand side
230           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
231           @param q: location of constraints           :param q: location of constraints
232           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
233           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
234           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
235           """           """
236           self.__domain=domain           self.__domain=domain
237           self.__D=D           self.__D=D
# Line 234  class NoPDE: Line 240  class NoPDE:
240           self.__r=r           self.__r=r
241           self.__u=None           self.__u=None
242           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
243    
244       def setReducedOn(self):       def setReducedOn(self):
245           """           """
246           sets the L{FunctionSpace} of the solution to L{ReducedSolution}           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
247           """           """
248           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
249           self.__u=None           self.__u=None
250    
251       def setReducedOff(self):       def setReducedOff(self):
252           """           """
253           sets the L{FunctionSpace} of the solution to L{Solution}           Sets the `FunctionSpace` of the solution to `Solution`.
254           """           """
255           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
256           self.__u=None           self.__u=None
257            
258       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
259           """           """
260           assigns values to the parameters.           Assigns values to the parameters.
261    
262           @param D: coefficient of the solution.           :param D: coefficient of the solution
263           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
264           @param Y: right hand side           :param Y: right hand side
265           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
266           @param q: location of constraints           :param q: location of constraints
267           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
268           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
269           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
270           """           """
271           if not D==None:           if not D==None:
272              self.__D=D              self.__D=D
# Line 276  class NoPDE: Line 283  class NoPDE:
283    
284       def getSolution(self):       def getSolution(self):
285           """           """
286           returns the solution           Returns the solution.
287            
288           @return: the solution of the problem           :return: the solution of the problem
289           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.           :rtype: `Data` object in the `FunctionSpace` `Solution` or
290                     `ReducedSolution`
291           """           """
292           if self.__u==None:           if self.__u==None:
293              if self.__D==None:              if self.__D==None:
# Line 296  class NoPDE: Line 304  class NoPDE:
304                  self.__u*=(1.-q)                  self.__u*=(1.-q)
305                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
306           return self.__u           return self.__u
307                
308  class Locator:  class Locator:
309       """       """
310       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
311       spatial coordinate x.         coordinate x.
312        
313       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
314       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.  
315       """       """
316    
317       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
318         """         """
319         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
320         or FunctionSpace where for the sample point which         or FunctionSpace for the sample point which is closest to the given
321         closest to the given point x.         point x.
322    
323         @param where: function space         :param where: function space
324         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
325         @param x: coefficient of the solution.         :param x: location(s) of the Locator
326         @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
327         """         """
328         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
329            self.__function_space=where            self.__function_space=where
330         else:         else:
331            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
332           iterative=False
333         if isinstance(x, list):         if isinstance(x, list):
334               if len(x)==0:
335                  raise "ValueError", "At least one point must be given."
336               try:
337                 iter(x[0])
338                 iterative=True
339               except TypeError:
340                 iterative=False
341           if iterative:
342             self.__id=[]             self.__id=[]
343             for p in x:             for p in x:
344                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 350  class Locator:
350         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
351         """         """
352         x=self.getX()         x=self.getX()
353         if instance(x,list):         if isinstance(x,list):
354            out="["            out="["
355            first=True            first=True
356            for xx in x:            for xx in x:
# Line 350  class Locator: Line 366  class Locator:
366    
367       def getX(self):       def getX(self):
368          """          """
369      Returns the exact coordinates of the Locator.          Returns the exact coordinates of the Locator.
370      """          """
371          return self(self.getFunctionSpace().getX())          return self(self.getFunctionSpace().getX())
372    
373       def getFunctionSpace(self):       def getFunctionSpace(self):
374          """          """
375      Returns the function space of the Locator.          Returns the function space of the Locator.
376      """          """
377          return self.__function_space          return self.__function_space
378    
379       def getId(self,item=None):       def getId(self,item=None):
380          """          """
381      Returns the identifier of the location.          Returns the identifier of the location.
382      """          """
383          if item == None:          if item == None:
384             return self.__id             return self.__id
385          else:          else:
# Line 375  class Locator: Line 391  class Locator:
391    
392       def __call__(self,data):       def __call__(self,data):
393          """          """
394      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.
395      the object is returned.          """
     """  
396          return self.getValue(data)          return self.getValue(data)
397    
398       def getValue(self,data):       def getValue(self,data):
399          """          """
400      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`
401      otherwise the object is returned.          object otherwise the object is returned.
402      """          """
403          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
404             if data.getFunctionSpace()==self.getFunctionSpace():             dat=util.interpolate(data,self.getFunctionSpace())
              dat=data  
            else:  
              dat=data.interpolate(self.getFunctionSpace())  
405             id=self.getId()             id=self.getId()
406             r=data.getRank()             r=data.getRank()
407             if isinstance(id,list):             if isinstance(id,list):
408                 out=[]                 out=[]
409                 for i in id:                 for i in id:
410                    o=data.getValueOfGlobalDataPoint(*i)                    o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
411                    if data.getRank()==0:                    if data.getRank()==0:
412                       out.append(o[0])                       out.append(o[0])
413                    else:                    else:
414                       out.append(o)                       out.append(o)
415                 return out                 return out
416             else:             else:
417               out=data.getValueOfGlobalDataPoint(*id)               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
418               if data.getRank()==0:               if data.getRank()==0:
419                  return out[0]                  return out[0]
420               else:               else:
# Line 412  class Locator: Line 424  class Locator:
424    
425  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
426     """     """
427     exceptions thrown by solvers     This is a generic exception thrown by solvers.
428     """     """
429     pass     pass
430    
431  class IndefinitePreconditioner(SolverSchemeException):  class IndefinitePreconditioner(SolverSchemeException):
432     """     """
433     the preconditioner is not positive definite.     Exception thrown if the preconditioner is not positive definite.
434     """     """
435     pass     pass
436    
437  class MaxIterReached(SolverSchemeException):  class MaxIterReached(SolverSchemeException):
438     """     """
439     maxium number of iteration steps is reached.     Exception thrown if the maximum number of iteration steps is reached.
440     """     """
441     pass     pass
442  class IterationBreakDown(SolverSchemeException):  
443    class CorrectionFailed(SolverSchemeException):
444     """     """
445     iteration scheme econouters an incurable breakdown.     Exception thrown if no convergence has been achieved in the solution
446       correction scheme.
447     """     """
448     pass     pass
449  class NegativeNorm(SolverSchemeException):  
450    class IterationBreakDown(SolverSchemeException):
451     """     """
452     a norm calculation returns a negative norm.     Exception thrown if the iteration scheme encountered an incurable breakdown.
453     """     """
454     pass     pass
455    
456  class IterationHistory(object):  class NegativeNorm(SolverSchemeException):
457     """     """
458     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.  
459     """     """
460     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]  
461    
462     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):
463         """     """
464         returns True if the C{norm_r} is C{tolerance}*C{norm_b}     Solver for
465    
466             *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}  
467    
468         """     with a symmetric and positive definite operator A (more details required!).
469         if TOL==None:     It uses the conjugate gradient method with preconditioner M providing an
470            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  
471    
472  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):     The iteration is terminated if
    """  
    Solver for  
473    
474     M{Ax=b}     *|r| <= atol+rtol*|r0|*
475    
476     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.  
477    
478     The iteration is terminated if the C{stoppingcriterium} function return C{True}.     *|r| = sqrt( bilinearform(Msolve(r),r))*
479    
480     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
481    
482     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,
483     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,
484     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst}.
485    
486     @param b: the right hand side of the liner system. C{b} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
487     @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)
488     @param Aprod: returns the value Ax     :param x: an initial guess for the solution
489     @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)
490     @param Msolve: solves Mx=r     :param Aprod: returns the value Ax
491     @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
492  type like argument C{x}.                  argument ``x``. The returned object needs to be of the same type
493     @param bilinearform: inner product C{<x,r>}                  like argument ``r``.
494     @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
495     @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
496     @type stoppingcriterium: function that returns C{True} or C{False}                   argument ``r``. The returned object needs to be of the same
497     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.                   type like argument ``x``.
498     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :param bilinearform: inner product ``<x,r>``
499     @param iter_max: maximum number of iteration steps.     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
500     @type iter_max: C{int}                         type like argument ``x`` and ``r`` is. The returned value
501     @return: the solution approximation and the corresponding residual                         is a ``float``.
502     @rtype: C{tuple}     :param atol: absolute tolerance
503     @warning: C{b} and C{x} are altered.     :type atol: non-negative ``float``
504       :param rtol: relative tolerance
505       :type rtol: non-negative ``float``
506       :param iter_max: maximum number of iteration steps
507       :type iter_max: ``int``
508       :return: the solution approximation and the corresponding residual
509       :rtype: ``tuple``
510       :warning: ``r`` and ``x`` are altered.
511     """     """
512     iter=0     iter=0
    if x==None:  
       x=0*b  
    else:  
       b += (-1)*Aprod(x)  
    r=b  
513     rhat=Msolve(r)     rhat=Msolve(r)
514     d = rhat     d = rhat
515     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
516     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm,"negative norm."
517       norm_r0=math.sqrt(rhat_dot_r)
518       atol2=atol+rtol*norm_r0
519       if atol2<=0:
520          raise ValueError,"Non-positive tolarance."
521       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
522    
523     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
524    
525    
526       while not math.sqrt(rhat_dot_r) <= atol2:
527         iter+=1         iter+=1
528         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
529    
530         q=Aprod(d)         q=Aprod(d)
531         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
532         x += alpha * d         x += alpha * d
533         r += (-alpha) * q         if isinstance(q,ArithmeticTuple):
534           r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
535           else:
536               r += (-alpha) * q
537         rhat=Msolve(r)         rhat=Msolve(r)
538         rhat_dot_r_new = bilinearform(rhat, r)         rhat_dot_r_new = bilinearform(rhat, r)
539         beta = rhat_dot_r_new / rhat_dot_r         beta = rhat_dot_r_new / rhat_dot_r
# Line 557  type like argument C{x}. Line 542  type like argument C{x}.
542    
543         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
544         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
545           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
546     return x,r     if verbose: print "PCG: tolerance reached after %s steps."%iter
547       return x,r,math.sqrt(rhat_dot_r)
548    
549  class Defect(object):  class Defect(object):
550      """      """
551      defines a non-linear defect F(x) of a variable x      Defines a non-linear defect F(x) of a variable x.
552      """      """
553      def __init__(self):      def __init__(self):
554          """          """
555          initialize defect          Initializes defect.
556          """          """
557          self.setDerivativeIncrementLength()          self.setDerivativeIncrementLength()
558    
559      def bilinearform(self, x0, x1):      def bilinearform(self, x0, x1):
560          """          """
561          returns the inner product of x0 and x1          Returns the inner product of x0 and x1
562          @param x0: a value for x  
563          @param x1: a value for x          :param x0: value for x0
564          @return: the inner product of x0 and x1          :param x1: value for x1
565          @rtype: C{float}          :return: the inner product of x0 and x1
566            :rtype: ``float``
567          """          """
568          return 0          return 0
569          
570      def norm(self,x):      def norm(self,x):
571          """          """
572          the norm of argument C{x}          Returns the norm of argument ``x``.
573    
574          @param x: a value for x          :param x: a value
575          @return: norm of argument x          :return: norm of argument x
576          @rtype: C{float}          :rtype: ``float``
577          @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
578          """          """
579          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
580          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
581          return math.sqrt(s)          return math.sqrt(s)
582    
   
583      def eval(self,x):      def eval(self,x):
584          """          """
585          returns the value F of a given x          Returns the value F of a given ``x``.
586    
587          @param x: value for which the defect C{F} is evalulated.          :param x: value for which the defect ``F`` is evaluated
588          @return: value of the defect at C{x}          :return: value of the defect at ``x``
589          """          """
590          return 0          return 0
591    
592      def __call__(self,x):      def __call__(self,x):
593          return self.eval(x)          return self.eval(x)
594    
595      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
596          """          """
597          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
598          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
599            direction of v with x as a starting point.
600    
601          @param inc: relative increment length          :param inc: relative increment length
602          @type inc: positive C{float}          :type inc: positive ``float``
603          """          """
604          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError,"positive increment required."
605          self.__inc=inc          self.__inc=inc
606    
607      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
608          """          """
609          returns the relative increment length used to approximate the derivative of the defect          Returns the relative increment length used to approximate the
610          @return: value of the defect at C{x}          derivative of the defect.
611          @rtype: positive C{float}          :return: value of the defect at ``x``
612            :rtype: positive ``float``
613          """          """
614          return self.__inc          return self.__inc
615    
616      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
617          """          """
618          returns the directional derivative at x0 in the direction of v          Returns the directional derivative at ``x0`` in the direction of ``v``.
619    
620          @param F0: value of this defect at x0          :param F0: value of this defect at x0
621          @param x0: value at which derivative is calculated.          :param x0: value at which derivative is calculated
622          @param v: direction          :param v: direction
623          @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
624          @return: derivative of this defect at x0 in the direction of C{v}                                  (self.norm(v)=0)
625          @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``
626          maybe oepsnew verwritten to use exact evalution.          :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
627                   used but this method maybe overwritten to use exact evaluation.
628          """          """
629          normx=self.norm(x0)          normx=self.norm(x0)
630          if normx>0:          if normx>0:
# Line 651  class Defect(object): Line 640  class Defect(object):
640          F1=self.eval(x0 + epsnew * v)          F1=self.eval(x0 + epsnew * v)
641          return (F1-F0)/epsnew          return (F1-F0)/epsnew
642    
643  ######################################      ######################################
644  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):
645     """     """
646     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
647       criterion:
648    
649     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
650      
651     where M{x0} is the initial guess.     where *x0* is the initial guess.
652    
653     @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
654     @type defect: L{Defect}                    *norm* used in the stopping criterion.
655     @param x: initial guess for the solution, C{x} is altered.     :type defect: `Defect`
656     @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.
657     @param iter_max: maximum number of iteration steps     :type x: any object type allowing basic operations such as
658     @type iter_max: positive C{int}              ``numpy.ndarray``, `Data`
659     @param sub_iter_max:     :param iter_max: maximum number of iteration steps
660     @type sub_iter_max:     :type iter_max: positive ``int``
661     @param atol: absolute tolerance for the solution     :param sub_iter_max: maximum number of inner iteration steps
662     @type atol: positive C{float}     :type sub_iter_max: positive ``int``
663     @param rtol: relative tolerance for the solution     :param atol: absolute tolerance for the solution
664     @type rtol: positive C{float}     :type atol: positive ``float``
665     @param gamma: tolerance safety factor for inner iteration     :param rtol: relative tolerance for the solution
666     @type gamma: positive C{float}, less than 1     :type rtol: positive ``float``
667     @param sub_tol_max: upper bound for inner tolerance.     :param gamma: tolerance safety factor for inner iteration
668     @type sub_tol_max: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
669     @return: an approximation of the solution with the desired accuracy     :param subtol_max: upper bound for inner tolerance
670     @rtype: same type as the initial guess.     :type subtol_max: positive ``float``, less than 1
671       :return: an approximation of the solution with the desired accuracy
672       :rtype: same type as the initial guess
673     """     """
674     lmaxit=iter_max     lmaxit=iter_max
675     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError,"atol needs to be non-negative."
676     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError,"rtol needs to be non-negative."
677     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."
678     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
679     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
680    
681     F=defect(x)     F=defect(x)
682     fnrm=defect.norm(F)     fnrm=defect.norm(F)
683     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
684     sub_tol=sub_tol_max     subtol=subtol_max
685     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
686     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print "             tolerance = %e."%subtol
687     iter=1     iter=1
688     #     #
689     # main iteration loop     # main iteration loop
690     #     #
691     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
692              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
693              #              #
694          #   adjust sub_tol_          #   adjust subtol_
695          #          #
696              if iter > 1:              if iter > 1:
697             rat=fnrm/fnrmo             rat=fnrm/fnrmo
698                 sub_tol_old=sub_tol                 subtol_old=subtol
699             sub_tol=gamma*rat**2             subtol=gamma*rat**2
700             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)
701             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
702          #          #
703          # calculate newton increment xc          # calculate newton increment xc
704              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
705              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
706              #     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
707              #              #
708              #                #
709              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
710              try:              try:
711                 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)
712              except MaxIterReached:              except MaxIterReached:
713                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
714              if sub_iter<0:              if sub_iter<0:
# Line 734  def NewtonGMRES(defect, x, iter_max=100, Line 726  def NewtonGMRES(defect, x, iter_max=100,
726    
727  def __givapp(c,s,vin):  def __givapp(c,s,vin):
728      """      """
729      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
730      @warning: C{vin} is altered.      ``vin``
731    
732        :warning: ``vin`` is altered.
733      """      """
734      vrot=vin      vrot=vin
735      if isinstance(c,float):      if isinstance(c,float):
736          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]]
737      else:      else:
738          for i in range(len(c)):          for i in range(len(c)):
739              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
740          w2=s[i]*vrot[i]+c[i]*vrot[i+1]          w2=s[i]*vrot[i]+c[i]*vrot[i+1]
741              vrot[i:i+2]=w1,w2              vrot[i]=w1
742                vrot[i+1]=w2
743      return vrot      return vrot
744    
745  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
746     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
747     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
748     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
749     g=numarray.zeros(iter_restart,numarray.Float64)     g=numpy.zeros(iter_restart,numpy.float64)
750     v=[]     v=[]
751    
752     rho=defect.norm(F0)     rho=defect.norm(F0)
753     if rho<=0.: return x0*0     if rho<=0.: return x0*0
754      
755     v.append(-F0/rho)     v.append(-F0/rho)
756     g[0]=rho     g[0]=rho
757     iter=0     iter=0
758     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
759            if iter  >= iter_max:
760      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached,"maximum number of %s steps reached."%iter_max
761    
762          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
763      v.append(p)          v.append(p)
764    
765      v_norm1=defect.norm(v[iter+1])          v_norm1=defect.norm(v[iter+1])
766    
767          # Modified Gram-Schmidt          # Modified Gram-Schmidt
768      for j in range(iter+1):          for j in range(iter+1):
769           h[j][iter]=defect.bilinearform(v[j],v[iter+1])                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
770           v[iter+1]-=h[j][iter]*v[j]              v[iter+1]-=h[j,iter]*v[j]
771          
772      h[iter+1][iter]=defect.norm(v[iter+1])          h[iter+1,iter]=defect.norm(v[iter+1])
773      v_norm2=h[iter+1][iter]          v_norm2=h[iter+1,iter]
774    
775          # Reorthogonalize if needed          # Reorthogonalize if needed
776      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)
777          for j in range(iter+1):                for j in range(iter+1):
778             hr=defect.bilinearform(v[j],v[iter+1])                  hr=defect.bilinearform(v[j],v[iter+1])
779                 h[j][iter]=h[j][iter]+hr                  h[j,iter]=h[j,iter]+hr
780                 v[iter+1] -= hr*v[j]                  v[iter+1] -= hr*v[j]
781    
782          v_norm2=defect.norm(v[iter+1])              v_norm2=defect.norm(v[iter+1])
783          h[iter+1][iter]=v_norm2              h[iter+1,iter]=v_norm2
784          #   watch out for happy breakdown          #   watch out for happy breakdown
785          if not v_norm2 == 0:          if not v_norm2 == 0:
786                  v[iter+1]=v[iter+1]/h[iter+1][iter]              v[iter+1]=v[iter+1]/h[iter+1,iter]
787    
788          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
789      if iter > 0 :          if iter > 0 :
790          hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numpy.zeros(iter+1,numpy.float64)
791          for i in range(iter+1) : hhat[i]=h[i][iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
792          hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
793              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
794    
795      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])
796    
797      if mu!=0 :          if mu!=0 :
798          c[iter]=h[iter][iter]/mu              c[iter]=h[iter,iter]/mu
799          s[iter]=-h[iter+1][iter]/mu              s[iter]=-h[iter+1,iter]/mu
800          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]
801          h[iter+1][iter]=0.0              h[iter+1,iter]=0.0
802          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])              gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
803                g[iter]=gg[0]
804                g[iter+1]=gg[1]
805    
806          # Update the residual norm          # Update the residual norm
807          rho=abs(g[iter+1])          rho=abs(g[iter+1])
808      iter+=1          iter+=1
809    
810     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
811     # It's time to compute x and leave.             # It's time to compute x and leave.
812     if iter > 0 :     if iter > 0 :
813       y=numarray.zeros(iter,numarray.Float64)           y=numpy.zeros(iter,numpy.float64)
814       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
815       if iter > 1 :         if iter > 1 :
816          i=iter-2            i=iter-2
817          while i>=0 :          while i>=0 :
818            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]
819            i=i-1            i=i-1
820       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
821       for i in range(iter-1):       for i in range(iter-1):
822      xhat += v[i]*y[i]      xhat += v[i]*y[i]
823     else :     else :
824        xhat=v[0] * 0        xhat=v[0] * 0
825    
826     if iter<iter_restart-1:     if iter<iter_restart-1:
827        stopped=iter        stopped=iter
828     else:     else:
829        stopped=-1        stopped=-1
830    
831     return xhat,stopped     return xhat,stopped
832    
833  ##############################################  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
834  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):     """
835  ################################################     Solver for
836    
837       *Ax=b*
838    
839       with a general operator A (more details required!).
840       It uses the generalized minimum residual method (GMRES).
841    
842       The iteration is terminated if
843    
844       *|r| <= atol+rtol*|r0|*
845    
846       where *r0* is the initial residual and *|.|* is the energy norm. In fact
847    
848       *|r| = sqrt( bilinearform(r,r))*
849    
850       :param r: initial residual *r=b-Ax*. ``r`` is altered.
851       :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
852       :param x: an initial guess for the solution
853       :type x: same like ``r``
854       :param Aprod: returns the value Ax
855       :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
856                    argument ``x``. The returned object needs to be of the same
857                    type like argument ``r``.
858       :param bilinearform: inner product ``<x,r>``
859       :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
860                           type like argument ``x`` and ``r``. The returned value is
861                           a ``float``.
862       :param atol: absolute tolerance
863       :type atol: non-negative ``float``
864       :param rtol: relative tolerance
865       :type rtol: non-negative ``float``
866       :param iter_max: maximum number of iteration steps
867       :type iter_max: ``int``
868       :param iter_restart: in order to save memory the orthogonalization process
869                            is terminated after ``iter_restart`` steps and the
870                            iteration is restarted.
871       :type iter_restart: ``int``
872       :return: the solution approximation and the corresponding residual
873       :rtype: ``tuple``
874       :warning: ``r`` and ``x`` are altered.
875       """
876     m=iter_restart     m=iter_restart
877       restarted=False
878     iter=0     iter=0
879     xc=x     if rtol>0:
880          r_dot_r = bilinearform(r, r)
881          if r_dot_r<0: raise NegativeNorm,"negative norm."
882          atol2=atol+rtol*math.sqrt(r_dot_r)
883          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
884       else:
885          atol2=atol
886          if verbose: print "GMRES: absolute tolerance = %e"%atol2
887       if atol2<=0:
888          raise ValueError,"Non-positive tolarance."
889    
890     while True:     while True:
891        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
892        xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)        if restarted:
893             r2 = r-Aprod(x-x2)
894          else:
895             r2=1*r
896          x2=x*1.
897          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
898          iter+=iter_restart
899        if stopped: break        if stopped: break
900        iter+=iter_restart            if verbose: print "GMRES: restart."
901     return xc        restarted=True
902       if verbose: print "GMRES: tolerance has been reached."
903       return x
904    
905  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):
906     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)  
907    
908     if x==None:     h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
909        x=0*b     c=numpy.zeros(iter_restart,numpy.float64)
910     else:     s=numpy.zeros(iter_restart,numpy.float64)
911        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)  
912     v=[]     v=[]
913    
914       r_dot_r = bilinearform(r, r)
915       if r_dot_r<0: raise NegativeNorm,"negative norm."
916     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
917      
918     v.append(r/rho)     v.append(r/rho)
919     g[0]=rho     g[0]=rho
920    
921     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
922       while not (rho<=atol or iter==iter_restart):
923    
924      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
925    
926      p=Msolve(Aprod(v[iter]))          if P_R!=None:
927                p=Aprod(P_R(v[iter]))
928            else:
929            p=Aprod(v[iter])
930      v.append(p)      v.append(p)
931    
932      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
933    
934    # Modified Gram-Schmidt
935        for j in range(iter+1):
936          h[j,iter]=bilinearform(v[j],v[iter+1])
937          v[iter+1]-=h[j,iter]*v[j]
938    
939  # Modified Gram-Schmidt      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
940      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]  
941    
942  # Reorthogonalize if needed  # Reorthogonalize if needed
943      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)
944       for j in range(iter+1):         for j in range(iter+1):
945          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
946              h[j][iter]=h[j][iter]+hr              h[j,iter]=h[j,iter]+hr
947              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
948    
949       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
950       h[iter+1][iter]=v_norm2       h[iter+1,iter]=v_norm2
951    
952  #   watch out for happy breakdown  #   watch out for happy breakdown
953          if not v_norm2 == 0:          if not v_norm2 == 0:
954           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
955    
956  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
957      if iter > 0 :      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
958          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])  
959    
960      if mu!=0 :      if mu!=0 :
961          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
962          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
963          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]
964          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
965          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
966                    g[iter]=gg[0]
967                    g[iter+1]=gg[1]
968  # Update the residual norm  # Update the residual norm
969                  
970          rho=abs(g[iter+1])          rho=abs(g[iter+1])
971            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
972      iter+=1      iter+=1
973    
974  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
975  # It's time to compute x and leave.          # It's time to compute x and leave.
976    
977     if iter > 0 :     if verbose: print "GMRES: iteration stopped after %s step."%iter
978       y=numarray.zeros(iter,numarray.Float64)         if iter > 0 :
979       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y=numpy.zeros(iter,numpy.float64)
980       if iter > 1 :         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
981          i=iter-2         if iter > 1 :
982            i=iter-2
983          while i>=0 :          while i>=0 :
984            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]
985            i=i-1            i=i-1
986       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
987       for i in range(iter-1):       for i in range(iter-1):
988      xhat += v[i]*y[i]      xhat += v[i]*y[i]
989     else : xhat=v[0]     else:
990         xhat=v[0] * 0
991     x += xhat     if P_R!=None:
992     if iter<iter_restart-1:        x += P_R(xhat)
993        stopped=True     else:
994     else:        x += xhat
995       if iter<iter_restart-1:
996          stopped=True
997       else:
998        stopped=False        stopped=False
999    
1000     return x,stopped     return x,stopped
1001    
1002  #################################################  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1003  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):      """
1004  #################################################      Solver for
1005      #  
1006      #  minres solves the system of linear equations Ax = b      *Ax=b*
1007      #  where A is a symmetric matrix (possibly indefinite or singular)  
1008      #  and b is a given vector.      with a symmetric and positive definite operator A (more details required!).
1009      #        It uses the minimum residual method (MINRES) with preconditioner M
1010      #  "A" may be a dense or sparse matrix (preferably sparse!)      providing an approximation of A.
1011      #  or the name of a function such that  
1012      #               y = A(x)      The iteration is terminated if
1013      #  returns the product y = Ax for any given vector x.  
1014      #      *|r| <= atol+rtol*|r0|*
1015      #  "M" defines a positive-definite preconditioner M = C C'.  
1016      #  "M" may be a dense or sparse matrix (preferably sparse!)      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1017      #  or the name of a function such that  
1018      #  solves the system My = x for any given vector x.      *|r| = sqrt( bilinearform(Msolve(r),r))*
1019      #  
1020      #      For details on the preconditioned conjugate gradient method see the book:
1021        
1022        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1023        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1024        C. Romine, and H. van der Vorst}.
1025    
1026        :param r: initial residual *r=b-Ax*. ``r`` is altered.
1027        :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1028        :param x: an initial guess for the solution
1029        :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1030        :param Aprod: returns the value Ax
1031        :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1032                     argument ``x``. The returned object needs to be of the same
1033                     type like argument ``r``.
1034        :param Msolve: solves Mx=r
1035        :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1036                      argument ``r``. The returned object needs to be of the same
1037                      type like argument ``x``.
1038        :param bilinearform: inner product ``<x,r>``
1039        :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1040                            type like argument ``x`` and ``r`` is. The returned value
1041                            is a ``float``.
1042        :param atol: absolute tolerance
1043        :type atol: non-negative ``float``
1044        :param rtol: relative tolerance
1045        :type rtol: non-negative ``float``
1046        :param iter_max: maximum number of iteration steps
1047        :type iter_max: ``int``
1048        :return: the solution approximation and the corresponding residual
1049        :rtype: ``tuple``
1050        :warning: ``r`` and ``x`` are altered.
1051        """
1052      #------------------------------------------------------------------      #------------------------------------------------------------------
1053      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
1054      # y  =  beta1 P' v1,  where  P = C**(-1).      # y  =  beta1 P' v1,  where  P = C**(-1).
1055      # v is really P' v1.      # v is really P' v1.
1056      #------------------------------------------------------------------      #------------------------------------------------------------------
1057      if x==None:      r1    = r
1058        x=0*b      y = Msolve(r)
1059      else:      beta1 = bilinearform(y,r)
       b += (-1)*Aprod(x)  
1060    
     r1    = b  
     y = Msolve(b)  
     beta1 = bilinearform(y,b)  
   
1061      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
1062    
1063      #  If b = 0 exactly, stop with x = 0.      #  If r = 0 exactly, stop with x
1064      if beta1==0: return x*0.      if beta1==0: return x
1065    
1066      if beta1> 0:      if beta1> 0: beta1  = math.sqrt(beta1)
       beta1  = math.sqrt(beta1)        
1067    
1068      #------------------------------------------------------------------      #------------------------------------------------------------------
1069      # Initialize quantities.      # Initialize quantities.
# Line 1008  def MINRES(b, Aprod, Msolve, bilinearfor Line 1083  def MINRES(b, Aprod, Msolve, bilinearfor
1083      ynorm2 = 0      ynorm2 = 0
1084      cs     = -1      cs     = -1
1085      sn     = 0      sn     = 0
1086      w      = b*0.      w      = r*0.
1087      w2     = b*0.      w2     = r*0.
1088      r2     = r1      r2     = r1
1089      eps    = 0.0001      eps    = 0.0001
1090    
1091      #---------------------------------------------------------------------      #---------------------------------------------------------------------
1092      # Main iteration loop.      # Main iteration loop.
1093      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1094      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1095    
1096      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
1097          iter    = iter  +  1          iter    = iter  +  1
# Line 1035  def MINRES(b, Aprod, Msolve, bilinearfor Line 1110  def MINRES(b, Aprod, Msolve, bilinearfor
1110          #-----------------------------------------------------------------          #-----------------------------------------------------------------
1111          s = 1/beta                 # Normalize previous vector (in y).          s = 1/beta                 # Normalize previous vector (in y).
1112          v = s*y                    # v = vk if P = I          v = s*y                    # v = vk if P = I
1113        
1114          y      = Aprod(v)          y      = Aprod(v)
1115        
1116          if iter >= 2:          if iter >= 2:
1117            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1118    
1119          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1120          y      += (- alfa/beta)*r2          y      += (- alfa/beta)*r2
1121          r1     = r2          r1     = r2
1122          r2     = y          r2     = y
1123          y = Msolve(r2)          y = Msolve(r2)
# Line 1052  def MINRES(b, Aprod, Msolve, bilinearfor Line 1127  def MINRES(b, Aprod, Msolve, bilinearfor
1127    
1128          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1129          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1130            
1131          if iter==1:                 # Initialize a few things.          if iter==1:                 # Initialize a few things.
1132            gmax   = abs( alfa )      # alpha1            gmax   = abs( alfa )      # alpha1
1133            gmin   = gmax             # alpha1            gmin   = gmax             # alpha1
# Line 1060  def MINRES(b, Aprod, Msolve, bilinearfor Line 1135  def MINRES(b, Aprod, Msolve, bilinearfor
1135          # Apply previous rotation Qk-1 to get          # Apply previous rotation Qk-1 to get
1136          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1137          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1138        
1139          oldeps = epsln          oldeps = epsln
1140          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1141          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 1145  def MINRES(b, Aprod, Msolve, bilinearfor
1145          # Compute the next plane rotation Qk          # Compute the next plane rotation Qk
1146    
1147          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1148          gamma  = max(gamma,eps)          gamma  = max(gamma,eps)
1149          cs     = gbar / gamma             # ck          cs     = gbar / gamma             # ck
1150          sn     = beta / gamma             # sk          sn     = beta / gamma             # sk
1151          phi    = cs * phibar              # phik          phi    = cs * phibar              # phik
# Line 1078  def MINRES(b, Aprod, Msolve, bilinearfor Line 1153  def MINRES(b, Aprod, Msolve, bilinearfor
1153    
1154          # Update  x.          # Update  x.
1155    
1156          denom = 1/gamma          denom = 1/gamma
1157          w1    = w2          w1    = w2
1158          w2    = w          w2    = w
1159          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1160          x     +=  phi*w          x     +=  phi*w
1161    
# Line 1095  def MINRES(b, Aprod, Msolve, bilinearfor Line 1170  def MINRES(b, Aprod, Msolve, bilinearfor
1170    
1171          # Estimate various norms and test for convergence.          # Estimate various norms and test for convergence.
1172    
1173          Anorm  = math.sqrt( tnorm2 )          Anorm  = math.sqrt( tnorm2 )
1174          ynorm  = math.sqrt( ynorm2 )          ynorm  = math.sqrt( ynorm2 )
1175    
1176          rnorm  = phibar          rnorm  = phibar
1177    
1178      return x      return x
1179    
1180  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):
1181      """
1182      Solver for
1183    
1184      *Ax=b*
1185    
1186  # TFQMR solver for linear systems    with a general operator A (more details required!).
1187  #    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
 #  
 # initialization  
 #  
   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  
1188    
1189    r=Msolve(r)    The iteration is terminated if
1190    
1191      *|r| <= atol+rtol*|r0|*
1192    
1193      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1194    
1195      *|r| = sqrt( bilinearform(r,r))*
1196    
1197      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1198      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1199      :param x: an initial guess for the solution
1200      :type x: same like ``r``
1201      :param Aprod: returns the value Ax
1202      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1203                   argument ``x``. The returned object needs to be of the same type
1204                   like argument ``r``.
1205      :param bilinearform: inner product ``<x,r>``
1206      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1207                          type like argument ``x`` and ``r``. The returned value is
1208                          a ``float``.
1209      :param atol: absolute tolerance
1210      :type atol: non-negative ``float``
1211      :param rtol: relative tolerance
1212      :type rtol: non-negative ``float``
1213      :param iter_max: maximum number of iteration steps
1214      :type iter_max: ``int``
1215      :rtype: ``tuple``
1216      :warning: ``r`` and ``x`` are altered.
1217      """
1218    u1=0    u1=0
1219    u2=0    u2=0
1220    y1=0    y1=0
1221    y2=0    y2=0
1222    
1223    w = r    w = r
1224    y1 = r    y1 = r
1225    iter = 0    iter = 0
1226    d = 0    d = 0
1227        v = Aprod(y1)
   v = Msolve(Aprod(y1))  
1228    u1 = v    u1 = v
1229      
1230    theta = 0.0;    theta = 0.0;
1231    eta = 0.0;    eta = 0.0;
1232    tau = math.sqrt(bilinearform(r,r))    rho=bilinearform(r,r)
1233    error = [ error, tau ]    if rho < 0: raise NegativeNorm,"negative norm."
1234    rho = tau * tau    tau = math.sqrt(rho)
1235    m=1    norm_r0=tau
1236  #    while tau>atol+rtol*norm_r0:
 #  TFQMR iteration  
 #  
 #  while ( iter < kmax-1 ):  
     
   while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):  
1237      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
1238    
1239      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1240        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
     if ( sigma == 0.0 ):  
       raise 'TFQMR breakdown, sigma=0'  
       
1241    
1242      alpha = rho / sigma      alpha = rho / sigma
1243    
# Line 1162  def TFQMR(b, Aprod, Msolve, bilinearform Line 1247  def TFQMR(b, Aprod, Msolve, bilinearform
1247  #  #
1248        if ( j == 1 ):        if ( j == 1 ):
1249          y2 = y1 - alpha * v          y2 = y1 - alpha * v
1250          u2 = Msolve(Aprod(y2))          u2 = Aprod(y2)
1251    
1252        m = 2 * (iter+1) - 2 + (j+1)        m = 2 * (iter+1) - 2 + (j+1)
1253        if j==0:        if j==0:
1254           w = w - alpha * u1           w = w - alpha * u1
1255           d = y1 + ( theta * theta * eta / alpha ) * d           d = y1 + ( theta * theta * eta / alpha ) * d
1256        if j==1:        if j==1:
# Line 1180  def TFQMR(b, Aprod, Msolve, bilinearform Line 1265  def TFQMR(b, Aprod, Msolve, bilinearform
1265  #  #
1266  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1267  #  #
1268       # 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'  
       
1269    
1270      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1271      beta = rhon / rho;      beta = rhon / rho;
1272      rho = rhon;      rho = rhon;
1273      y1 = w + beta * y2;      y1 = w + beta * y2;
1274      u1 = Msolve(Aprod(y1))      u1 = Aprod(y1)
1275      v = u1 + beta * ( u2 + beta * v )      v = u1 + beta * ( u2 + beta * v )
1276      error = [ error, tau ]  
1277      total_iters = iter      iter += 1
       
     iter = iter + 1  
1278    
1279    return x    return x
1280    
# Line 1208  def TFQMR(b, Aprod, Msolve, bilinearform Line 1283  def TFQMR(b, Aprod, Msolve, bilinearform
1283    
1284  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1285     """     """
1286     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
1287       ArithmeticTuple and ``a`` is a float.
1288    
1289     example of usage:     Example of usage::
1290    
1291     from esys.escript import Data         from esys.escript import Data
1292     from numarray import array         from numpy import array
1293     a=Data(...)         a=Data(...)
1294     b=array([1.,4.])         b=array([1.,4.])
1295     x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1296     y=5.*x         y=5.*x
1297    
1298     """     """
1299     def __init__(self,*args):     def __init__(self,*args):
1300         """         """
1301         initialize object with elements args.         Initializes object with elements ``args``.
1302    
1303         @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
1304                        scaling (x=a*y)
1305         """         """
1306         self.__items=list(args)         self.__items=list(args)
1307    
1308     def __len__(self):     def __len__(self):
1309         """         """
1310         number of items         Returns the number of items.
1311    
1312         @return: number of items         :return: number of items
1313         @rtype: C{int}         :rtype: ``int``
1314         """         """
1315         return len(self.__items)         return len(self.__items)
1316    
1317     def __getitem__(self,index):     def __getitem__(self,index):
1318         """         """
1319         get an item         Returns item at specified position.
1320    
1321         @param index: item to be returned         :param index: index of item to be returned
1322         @type index: C{int}         :type index: ``int``
1323         @return: item with index C{index}         :return: item with index ``index``
1324         """         """
1325         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1326    
1327     def __mul__(self,other):     def __mul__(self,other):
1328         """         """
1329         scaling from the right         Scales by ``other`` from the right.
1330    
1331         @param other: scaling factor         :param other: scaling factor
1332         @type other: C{float}         :type other: ``float``
1333         @return: itemwise self*other         :return: itemwise self*other
1334         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1335         """         """
1336         out=[]         out=[]
1337         try:           try:
1338             l=len(other)             l=len(other)
1339             if l!=len(self):             if l!=len(self):
1340                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1341             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1342         except TypeError:         except TypeError:
1343         for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
1344         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1345    
1346     def __rmul__(self,other):     def __rmul__(self,other):
1347         """         """
1348         scaling from the left         Scales by ``other`` from the left.
1349    
1350         @param other: scaling factor         :param other: scaling factor
1351         @type other: C{float}         :type other: ``float``
1352         @return: itemwise other*self         :return: itemwise other*self
1353         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1354         """         """
1355         out=[]         out=[]
1356         try:           try:
1357             l=len(other)             l=len(other)
1358             if l!=len(self):             if l!=len(self):
1359                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1360             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1361         except TypeError:         except TypeError:
1362         for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
1363         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1364    
1365     def __div__(self,other):     def __div__(self,other):
1366         """         """
1367         dividing from the right         Scales by (1/``other``) from the right.
1368    
1369         @param other: scaling factor         :param other: scaling factor
1370         @type other: C{float}         :type other: ``float``
1371         @return: itemwise self/other         :return: itemwise self/other
1372         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1373         """         """
1374         return self*(1/other)         return self*(1/other)
1375    
1376     def __rdiv__(self,other):     def __rdiv__(self,other):
1377         """         """
1378         dividing from the left         Scales by (1/``other``) from the left.
1379    
1380         @param other: scaling factor         :param other: scaling factor
1381         @type other: C{float}         :type other: ``float``
1382         @return: itemwise other/self         :return: itemwise other/self
1383         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1384         """         """
1385         out=[]         out=[]
1386         try:           try:
1387             l=len(other)             l=len(other)
1388             if l!=len(self):             if l!=len(self):
1389                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1390             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1391         except TypeError:         except TypeError:
1392         for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
1393         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1394      
1395     def __iadd__(self,other):     def __iadd__(self,other):
1396         """         """
1397         in-place add of other to self         Inplace addition of ``other`` to self.
1398    
1399         @param other: increment         :param other: increment
1400         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1401         """         """
1402         if len(self) != len(other):         if len(self) != len(other):
1403             raise ValueError,"tuple length must match."             raise ValueError,"tuple lengths must match."
1404         for i in range(len(self)):         for i in range(len(self)):
1405             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1406         return self         return self
1407    
1408     def __add__(self,other):     def __add__(self,other):
1409         """         """
1410         add other to self         Adds ``other`` to self.
1411    
1412         @param other: increment         :param other: increment
1413         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1414         """         """
1415         out=[]         out=[]
1416         try:           try:
1417             l=len(other)             l=len(other)
1418             if l!=len(self):             if l!=len(self):
1419                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1420             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1421         except TypeError:         except TypeError:
1422         for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
1423         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1424    
1425     def __sub__(self,other):     def __sub__(self,other):
1426         """         """
1427         subtract other from self         Subtracts ``other`` from self.
1428    
1429         @param other: increment         :param other: decrement
1430         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1431         """         """
1432         out=[]         out=[]
1433         try:           try:
1434             l=len(other)             l=len(other)
1435             if l!=len(self):             if l!=len(self):
1436                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1437             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1438         except TypeError:         except TypeError:
1439         for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
1440         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1441      
1442     def __isub__(self,other):     def __isub__(self,other):
1443         """         """
1444         subtract other from self         Inplace subtraction of ``other`` from self.
1445    
1446         @param other: increment         :param other: decrement
1447         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1448         """         """
1449         if len(self) != len(other):         if len(self) != len(other):
1450             raise ValueError,"tuple length must match."             raise ValueError,"tuple length must match."
# Line 1377  class ArithmeticTuple(object): Line 1454  class ArithmeticTuple(object):
1454    
1455     def __neg__(self):     def __neg__(self):
1456         """         """
1457         negate         Negates values.
   
1458         """         """
1459         out=[]         out=[]
1460         for i in range(len(self)):         for i in range(len(self)):
# Line 1388  class ArithmeticTuple(object): Line 1464  class ArithmeticTuple(object):
1464    
1465  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1466        """        """
1467        This provides a framwork for solving linear homogeneous saddle point problem of the form        This class provides a framework for solving linear homogeneous saddle
1468          point problems of the form::
              Av+B^*p=f  
              Bv    =0  
1469    
1470        for the unknowns v and p and given operators A and B and given right hand side f.            *Av+B^*p=f*
1471        B^* is the adjoint operator of B is the given inner product.            *Bv     =0*
1472    
1473          for the unknowns *v* and *p* and given operators *A* and *B* and
1474          given right hand side *f*. *B^** is the adjoint operator of *B*.
1475          *A* may depend weakly on *v* and *p*.
1476        """        """
1477        def __init__(self,**kwargs):        def __init__(self, **kwargs):
1478        """
1479        initializes the saddle point problem
1480        """
1481            self.resetControlParameters()
1482          self.setTolerance()          self.setTolerance()
1483          self.setToleranceReductionFactor()          self.setAbsoluteTolerance()
1484          def resetControlParameters(self,gamma=0.85, gamma_min=1.e-8,chi_max=0.1, omega_div=0.2, omega_conv=1.1, rtol_min=1.e-7, rtol_max=0.1, chi=1., C_p=1., C_v=1., safety_factor=0.3):
1485             """
1486             sets a control parameter
1487    
1488             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1489             :type gamma: ``float``
1490             :param gamma_min: minimum value for ``gamma``.
1491             :type gamma_min: ``float``
1492             :param chi_max: maximum tolerable converegence rate.
1493             :type chi_max: ``float``
1494             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1495             :type omega_div: ``float``
1496             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1497             :type omega_conv: ``float``
1498             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1499             :type rtol_min: ``float``
1500             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1501             :type rtol_max: ``float``
1502             :param chi: initial convergence measure.
1503             :type chi: ``float``
1504             :param C_p: initial value for constant to adjust pressure tolerance
1505             :type C_p: ``float``
1506             :param C_v: initial value for constant to adjust velocity tolerance
1507             :type C_v: ``float``
1508             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1509             :type safety_factor: ``float``
1510             """
1511             self.setControlParameter(gamma, gamma_min ,chi_max , omega_div , omega_conv, rtol_min , rtol_max, chi,C_p, C_v,safety_factor)
1512    
1513          def setControlParameter(self,gamma=None, gamma_min=None ,chi_max=None, omega_div=None, omega_conv=None, rtol_min=None, rtol_max=None, chi=None, C_p=None, C_v=None, safety_factor=None):
1514             """
1515             sets a control parameter
1516    
1517             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1518             :type gamma: ``float``
1519             :param gamma_min: minimum value for ``gamma``.
1520             :type gamma_min: ``float``
1521             :param chi_max: maximum tolerable converegence rate.
1522             :type chi_max: ``float``
1523             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1524             :type omega_div: ``float``
1525             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1526             :type omega_conv: ``float``
1527             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1528             :type rtol_min: ``float``
1529             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1530             :type rtol_max: ``float``
1531             :param chi: initial convergence measure.
1532             :type chi: ``float``
1533             :param C_p: initial value for constant to adjust pressure tolerance
1534             :type C_p: ``float``
1535             :param C_v: initial value for constant to adjust velocity tolerance
1536             :type C_v: ``float``
1537             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1538             :type safety_factor: ``float``
1539             """
1540             if not gamma == None:
1541                if gamma<=0 or gamma>=1:
1542                   raise ValueError,"Initial gamma needs to be positive and less than 1."
1543             else:
1544                gamma=self.__gamma
1545    
1546             if not gamma_min == None:
1547                if gamma_min<=0 or gamma_min>=1:
1548                   raise ValueError,"gamma_min needs to be positive and less than 1."
1549             else:
1550                gamma_min = self.__gamma_min
1551    
1552             if not chi_max == None:
1553                if chi_max<=0 or chi_max>=1:
1554                   raise ValueError,"chi_max needs to be positive and less than 1."
1555             else:
1556                chi_max = self.__chi_max
1557    
1558             if not omega_div == None:
1559                if omega_div<=0 or omega_div >=1:
1560                   raise ValueError,"omega_div needs to be positive and less than 1."
1561             else:
1562                omega_div=self.__omega_div
1563    
1564             if not omega_conv == None:
1565                if omega_conv<1:
1566                   raise ValueError,"omega_conv needs to be greater or equal to 1."
1567             else:
1568                omega_conv=self.__omega_conv
1569    
1570             if not rtol_min == None:
1571                if rtol_min<=0 or rtol_min>=1:
1572                   raise ValueError,"rtol_min needs to be positive and less than 1."
1573             else:
1574                rtol_min=self.__rtol_min
1575    
1576             if not rtol_max == None:
1577                if rtol_max<=0 or rtol_max>=1:
1578                   raise ValueError,"rtol_max needs to be positive and less than 1."
1579             else:
1580                rtol_max=self.__rtol_max
1581    
1582             if not chi == None:
1583                if chi<=0 or chi>1:
1584                   raise ValueError,"chi needs to be positive and less than 1."
1585             else:
1586                chi=self.__chi
1587    
1588             if not C_p == None:
1589                if C_p<1:
1590                   raise ValueError,"C_p need to be greater or equal to 1."
1591             else:
1592                C_p=self.__C_p
1593    
1594             if not C_v == None:
1595                if C_v<1:
1596                   raise ValueError,"C_v need to be greater or equal to 1."
1597             else:
1598                C_v=self.__C_v
1599    
1600             if not safety_factor == None:
1601                if safety_factor<=0 or safety_factor>1:
1602                   raise ValueError,"safety_factor need to be between zero and one."
1603             else:
1604                safety_factor=self.__safety_factor
1605    
1606             if gamma<gamma_min:
1607                   raise ValueError,"gamma = %e needs to be greater or equal gamma_min = %e."%(gamma,gamma_min)
1608             if rtol_max<=rtol_min:
1609                   raise ValueError,"rtol_max = %e needs to be greater rtol_min = %e."%(rtol_max,rtol_min)
1610                
1611             self.__gamma = gamma
1612             self.__gamma_min = gamma_min
1613             self.__chi_max = chi_max
1614             self.__omega_div = omega_div
1615             self.__omega_conv = omega_conv
1616             self.__rtol_min = rtol_min
1617             self.__rtol_max = rtol_max
1618             self.__chi = chi
1619             self.__C_p = C_p
1620             self.__C_v = C_v
1621             self.__safety_factor = safety_factor
1622    
1623          #=============================================================
1624        def initialize(self):        def initialize(self):
1625          """          """
1626          initialize the problem (overwrite)          Initializes the problem (overwrite).
1627          """          """
1628          pass          pass
       def B(self,v):  
          """  
          returns Bv (overwrite)  
          @rtype: equal to the type of p  
1629    
1630           @note: boundary conditions on p should be zero!        def inner_pBv(self,p,Bv):
1631           """           """
1632           pass           Returns inner product of element p and Bv (overwrite).
1633    
1634        def inner(self,p0,p1):           :param p: a pressure increment
1635             :param Bv: a residual
1636             :return: inner product of element p and Bv
1637             :rtype: ``float``
1638             :note: used if PCG is applied.
1639           """           """
1640           returns inner product of two element p0 and p1  (overwrite)           raise NotImplementedError,"no inner product for p and Bv implemented."
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
1641    
1642           @rtype: equal to the type of p        def inner_p(self,p0,p1):
1643           """           """
1644           pass           Returns inner product of p0 and p1 (overwrite).
1645    
1646        def solve_A(self,u,p):           :param p0: a pressure
1647             :param p1: a pressure
1648             :return: inner product of p0 and p1
1649             :rtype: ``float``
1650             """
1651             raise NotImplementedError,"no inner product for p implemented."
1652      
1653          def norm_v(self,v):
1654           """           """
1655           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           Returns the norm of v (overwrite).
1656    
1657           @rtype: equal to the type of v           :param v: a velovity
1658           @note: boundary conditions on v should be zero!           :return: norm of v
1659             :rtype: non-negative ``float``
1660             """
1661             raise NotImplementedError,"no norm of v implemented."
1662          def getDV(self, p, v, tol):
1663           """           """
1664           pass           return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
1665    
1666        def solve_prec(self,p):           :param p: pressure
1667             :param v: pressure
1668             :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1669             :note: Only *A* may depend on *v* and *p*
1670           """           """
1671           provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)           raise NotImplementedError,"no dv calculation implemented."
1672    
1673            
1674          def Bv(self,v, tol):
1675            """
1676            Returns Bv with accuracy `tol` (overwrite)
1677    
1678            :rtype: equal to the type of p
1679            :note: boundary conditions on p should be zero!
1680            """
1681            raise NotImplementedError, "no operator B implemented."
1682    
1683           @rtype: equal to the type of p        def norm_Bv(self,Bv):
1684            """
1685            Returns the norm of Bv (overwrite).
1686    
1687            :rtype: equal to the type of p
1688            :note: boundary conditions on p should be zero!
1689            """
1690            raise NotImplementedError, "no norm of Bv implemented."
1691    
1692          def solve_AinvBt(self,dp, tol):
1693           """           """
1694           pass           Solves *A dv=B^*dp* with accuracy `tol`
1695    
1696        def stoppingcriterium(self,Bv,v,p):           :param dp: a pressure increment
1697             :return: the solution of *A dv=B^*dp*
1698             :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1699           """           """
1700           returns a True if iteration is terminated. (overwrite)           raise NotImplementedError,"no operator A implemented."
1701    
1702           @rtype: C{bool}        def solve_prec(self,Bv, tol):
1703           """           """
1704           pass           Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
               
       def __inner(self,p,r):  
          return self.inner(p,r[1])  
1705    
1706        def __inner_p(self,p1,p2):           :rtype: equal to the type of p
1707           return self.inner(p1,p2)           :note: boundary conditions on p should be zero!
1708                   """
1709        def __inner_a(self,a1,a2):           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1710           return self.inner_a(a1,a2)        #=============================================================
1711          def __Aprod_PCG(self,dp):
1712              dv=self.solve_AinvBt(dp, self.__subtol)
1713              return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1714    
1715          def __inner_PCG(self,p,r):
1716             return self.inner_pBv(p,r[1])
1717    
1718          def __Msolve_PCG(self,r):
1719              return self.solve_prec(r[1], self.__subtol)
1720          #=============================================================
1721          def __Aprod_GMRES(self,p):
1722              return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1723          def __inner_GMRES(self,p0,p1):
1724             return self.inner_p(p0,p1)
1725    
1726          #=============================================================
1727          def norm_p(self,p):
1728              """
1729              calculates the norm of ``p``
1730              
1731              :param p: a pressure
1732              :return: the norm of ``p`` using the inner product for pressure
1733              :rtype: ``float``
1734              """
1735              f=self.inner_p(p,p)
1736              if f<0: raise ValueError,"negative pressure norm."
1737              return math.sqrt(f)
1738          
1739          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1740             """
1741             Solves the saddle point problem using initial guesses v and p.
1742    
1743        def __inner_a1(self,a1,a2):           :param v: initial guess for velocity
1744           return self.inner(a1[1],a2[1])           :param p: initial guess for pressure
1745             :type v: `Data`
1746             :type p: `Data`
1747             :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1748             :param max_iter: maximum number of iteration steps per correction
1749                              attempt
1750             :param verbose: if True, shows information on the progress of the
1751                             saddlepoint problem solver.
1752             :param iter_restart: restart the iteration after ``iter_restart`` steps
1753                                  (only used if useUzaw=False)
1754             :type usePCG: ``bool``
1755             :type max_iter: ``int``
1756             :type verbose: ``bool``
1757             :type iter_restart: ``int``
1758             :rtype: ``tuple`` of `Data` objects
1759             """
1760             self.verbose=verbose
1761             rtol=self.getTolerance()
1762             atol=self.getAbsoluteTolerance()
1763             correction_step=0
1764             converged=False
1765             error=None
1766             chi=None
1767             gamma=self.__gamma
1768             C_p=self.__C_p
1769             C_v=self.__C_v
1770             while not converged:
1771                  if error== None or chi == None:
1772                      gamma_new=gamma/self.__omega_conv
1773                  else:
1774                     if chi < self.__chi_max:
1775                        gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)
1776                     else:
1777                        gamma_new=gamma*self.__omega_div
1778                  gamma=max(gamma_new, self.__gamma_min)
1779                  # calculate velocity for current pressure:
1780                  rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)
1781                  rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)
1782                  self.__subtol=rtol_p**2
1783                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: gamma = %e, rtol_v= %e, rtol_p=%e"%(correction_step,gamma,rtol_v,rtol_p)
1784                  if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol
1785                  # calculate velocity for current pressure: A*dv= F-A*v-B*p
1786                  dv1=self.getDV(p,v,rtol_v)
1787                  v1=v+dv1
1788                  Bv1=self.Bv(v1, self.__subtol)
1789                  norm_Bv1=self.norm_Bv(Bv1)
1790                  norm_dv1=self.norm_v(dv1)
1791                  norm_v1=self.norm_v(v1)
1792                  ATOL=norm_v1*rtol+atol
1793                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv1, norm_dv1, norm_v1)
1794                  if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1795                  if max(norm_Bv1, norm_dv1) <= ATOL:
1796                      converged=True
1797                      v=v1
1798                  else:
1799                      # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1800                      if usePCG:
1801                        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)
1802                        v2=r[0]
1803                        Bv2=r[1]
1804                      else:
1805                        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)
1806                        dv2=self.solve_AinvBt(dp, self.__subtol)
1807                        v2=v1-dv2
1808                        Bv2=self.Bv(v2, self.__subtol)
1809                      #
1810                      # convergence indicators:
1811                      #
1812                      norm_v2=self.norm_v(v2)
1813                      norm_dv2=self.norm_v(v2-v)
1814                      error_new=max(norm_dv2, norm_Bv1)
1815                      norm_Bv2=self.norm_Bv(Bv2)
1816                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s (part 2): Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv2, norm_dv2, norm_v2)
1817                      if error !=None:
1818                          chi_new=error_new/error
1819                          if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e"%(correction_step,chi_new)
1820                          if chi != None:
1821                              gamma0=max(gamma, 1-chi/chi_new)
1822                              C_p*=gamma0/gamma
1823                              C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)
1824                          chi=chi_new
1825                      error = error_new
1826                      correction_step+=1
1827                      if correction_step>max_correction_steps:
1828                            raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1829                      v,p=v2,p+dp
1830             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1831         return v,p
1832    
1833        def __stoppingcriterium(self,norm_r,r,p):        #========================================================================
1834            return self.stoppingcriterium(r[1],r[0],p)        def setTolerance(self,tolerance=1.e-4):
1835             """
1836             Sets the relative tolerance for (v,p).
1837    
1838        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):           :param tolerance: tolerance to be used
1839            return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)           :type tolerance: non-negative ``float``
1840             """
1841             if tolerance<0:
1842                 raise ValueError,"tolerance must be positive."
1843             self.__rtol=tolerance
1844    
       def setTolerance(self,tolerance=1.e-8):  
               self.__tol=tolerance  
1845        def getTolerance(self):        def getTolerance(self):
1846                return self.__tol           """
1847        def setToleranceReductionFactor(self,reduction=0.01):           Returns the relative tolerance.
1848                self.__reduction=reduction  
1849        def getSubProblemTolerance(self):           :return: relative tolerance
1850                return self.__reduction*self.getTolerance()           :rtype: ``float``
1851             """
1852        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):           return self.__rtol
1853                """  
1854                solves the saddle point problem using initial guesses v and p.        def setAbsoluteTolerance(self,tolerance=0.):
1855             """
1856                @param max_iter: maximum number of iteration steps.           Sets the absolute tolerance.
1857                """  
1858                self.verbose=verbose           :param tolerance: tolerance to be used
1859                self.show_details=show_details and self.verbose           :type tolerance: non-negative ``float``
1860             """
1861                # assume p is known: then v=A^-1(f-B^*p)           if tolerance<0:
1862                # which leads to BA^-1B^*p = BA^-1f                 raise ValueError,"tolerance must be non-negative."
1863             self.__atol=tolerance
1864            # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)        
1865            self.__z=v+self.solve_A(v,p*0)        def getAbsoluteTolerance(self):
1866                Bz=self.B(self.__z)           """
1867                #           Returns the absolute tolerance.
1868            #   solve BA^-1B^*p = Bz  
1869                #           :return: absolute tolerance
1870                #           :rtype: ``float``
1871                #           """
1872                self.iter=0           return self.__atol
           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)  
1873    
1874    
1875  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1876     """     """
1877     creates a mask on the Solution(domain) function space which one for samples     Creates a mask on the Solution(domain) function space where the value is
1878     that touch the boundary tagged by tags.     one for samples that touch the boundary tagged by tags.
1879    
1880     usage: m=MaskFromBoundaryTag(domain,"left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1881    
1882     @param domain: a given domain     :param domain: domain to be used
1883     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1884     @param tags: boundray tags     :param tags: boundary tags
1885     @type tags: C{str}     :type tags: ``str``
1886     @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.     :return: a mask which marks samples that are touching the boundary tagged
1887     @rtype: L{escript.Data} of rank 0              by any of the given tags
1888       :rtype: `escript.Data` of rank 0
1889     """     """
1890     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1891     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1892     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1893     pde.setValue(y=d)     pde.setValue(y=d)
1894     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
 #============================================================================================================================================  
 class SaddlePointProblem(object):  
    """  
    This implements a solver for a saddlepoint problem  
   
    M{f(u,p)=0}  
    M{g(u)=0}  
   
    for u and p. The problem is solved with an inexact Uszawa scheme for p:  
1895    
1896     M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}  def MaskFromTag(domain,*tags):
    M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}  
   
    where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of  
    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'  
    Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays  
    in fact the role of a preconditioner.  
1897     """     """
1898     def __init__(self,verbose=False,*args):     Creates a mask on the Solution(domain) function space where the value is
1899         """     one for samples that touch regions tagged by tags.
        initializes the problem  
1900    
1901         @param verbose: switches on the printing out some information     Usage: m=MaskFromTag(domain, "ham")
        @type verbose: C{bool}  
        @note: this method may be overwritten by a particular saddle point problem  
        """  
        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.")  
   
    def trace(self,text):  
        """  
        prints text if verbose has been set  
   
        @param text: a text message  
        @type text: C{str}  
        """  
        if self.__verbose: print "%s: %s"%(str(self),text)  
   
    def solve_f(self,u,p,tol=1.e-8):  
        """  
        solves  
1902    
1903         A_f du = f(u,p)     :param domain: domain to be used
1904       :type domain: `escript.Domain`
1905         with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.     :param tags: boundary tags
1906       :type tags: ``str``
1907         @param u: current approximation of u     :return: a mask which marks samples that are touching the boundary tagged
1908         @type u: L{escript.Data}              by any of the given tags
1909         @param p: current approximation of p     :rtype: `escript.Data` of rank 0
1910         @type p: L{escript.Data}     """
1911         @param tol: tolerance expected for du     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1912         @type tol: C{float}     d=escript.Scalar(0.,escript.Function(domain))
1913         @return: increment du     for t in tags: d.setTaggedValue(t,1.)
1914         @rtype: L{escript.Data}     pde.setValue(Y=d)
1915         @note: this method has to be overwritten by a particular saddle point problem     return util.whereNonZero(pde.getRightHandSide())
        """  
        pass  
   
    def solve_g(self,u,tol=1.e-8):  
        """  
        solves  
   
        Q_g dp = g(u)  
   
        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.  
   
        @param u: current approximation of u  
        @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  
   
    def inner(self,p0,p1):  
        """  
        inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)  
        @return: inner product of p0 and p1  
        @rtype: C{float}  
        """  
        pass  
1916    
    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  
1917    
         @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.2719

  ViewVC Help
Powered by ViewVC 1.1.26