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

  ViewVC Help
Powered by ViewVC 1.1.26