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

  ViewVC Help
Powered by ViewVC 1.1.26