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

  ViewVC Help
Powered by ViewVC 1.1.26