/[escript]/trunk/escript/py_src/util.py
ViewVC logotype

Diff of /trunk/escript/py_src/util.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 775 by ksteube, Mon Jul 10 04:00:08 2006 UTC revision 1042 by gross, Mon Mar 19 03:50:34 2007 UTC
# Line 26  import math Line 26  import math
26  import numarray  import numarray
27  import escript  import escript
28  import os  import os
29    from esys.escript import C_GeneralTensorProduct
30    
31  #=========================================================  #=========================================================
32  #   some helpers:  #   some helpers:
33  #=========================================================  #=========================================================
34    def getTagNames(domain):
35        """
36        returns a list of the tag names used by the domain
37    
38        
39        @param domain: a domain object
40        @type domain: C{escript.Domain}
41        @return: a list of the tag name used by the domain.
42        @rtype: C{list} of C{str}
43        """
44        return [n.strip() for n in domain.showTagNames().split(",") ]
45        
46  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
47      """      """
48      writes a L{Data} objects into a files using the the VTK XML file format.      writes a L{Data} objects into a files using the the VTK XML file format.
# Line 238  def inf(arg): Line 251  def inf(arg):
251  #=========================================================================  #=========================================================================
252  #   some little helpers  #   some little helpers
253  #=========================================================================  #=========================================================================
254  def pokeShape(arg):  def getRank(arg):
255        """
256        identifies the rank of its argument
257    
258        @param arg: a given object
259        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
260        @return: the rank of the argument
261        @rtype: C{int}
262        @raise TypeError: if type of arg cannot be processed
263        """
264    
265        if isinstance(arg,numarray.NumArray):
266            return arg.rank
267        elif isinstance(arg,escript.Data):
268            return arg.getRank()
269        elif isinstance(arg,float):
270            return 0
271        elif isinstance(arg,int):
272            return 0
273        elif isinstance(arg,Symbol):
274            return arg.getRank()
275        else:
276          raise TypeError,"getShape: cannot identify shape"
277    def getShape(arg):
278      """      """
279      identifies the shape of its argument      identifies the shape of its argument
280    
# Line 260  def pokeShape(arg): Line 296  def pokeShape(arg):
296      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
297          return arg.getShape()          return arg.getShape()
298      else:      else:
299        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
300    
301  def pokeDim(arg):  def pokeDim(arg):
302      """      """
# Line 283  def commonShape(arg0,arg1): Line 319  def commonShape(arg0,arg1):
319      """      """
320      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.
321    
322      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
323      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
324      @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.      @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.
325      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
326      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
327      """      """
328      sh0=pokeShape(arg0)      sh0=getShape(arg0)
329      sh1=pokeShape(arg1)      sh1=getShape(arg1)
330      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
331         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
332               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
# Line 440  def matchShape(arg0,arg1): Line 476  def matchShape(arg0,arg1):
476      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
477      @param arg1: a given object      @param arg1: a given object
478      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
479      @return: arg0 and arg1 where copies are returned when the shape has to be changed.      @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
480      @rtype: C{tuple}      @rtype: C{tuple}
481      """      """
482      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
483      sh0=pokeShape(arg0)      sh0=getShape(arg0)
484      sh1=pokeShape(arg1)      sh1=getShape(arg1)
485      if len(sh0)<len(sh):      if len(sh0)<len(sh):
486         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
487      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
# Line 573  class Symbol(object): Line 609  class Symbol(object):
609            if isinstance(a,Symbol):            if isinstance(a,Symbol):
610               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
611            else:            else:
612                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
613                if len(s)>0:                if len(s)>0:
614                   out.append(numarray.zeros(s),numarray.Float64)                   out.append(numarray.zeros(s),numarray.Float64)
615                else:                else:
# Line 1310  def whereNonZero(arg,tol=0.): Line 1346  def whereNonZero(arg,tol=0.):
1346     else:     else:
1347        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1348    
1349    def erf(arg):
1350       """
1351       returns erf of argument arg
1352    
1353       @param arg: argument
1354       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1355       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1356       @raises TypeError: if the type of the argument is not expected.
1357       """
1358       if isinstance(arg,escript.Data):
1359          return arg._erf()
1360       else:
1361          raise TypeError,"erf: Unknown argument type."
1362    
1363  def sin(arg):  def sin(arg):
1364     """     """
1365     returns sine of argument arg     returns sine of argument arg
# Line 2930  def trace(arg,axis_offset=0): Line 2980  def trace(arg,axis_offset=0):
2980    
2981     @param arg: argument     @param arg: argument
2982     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2983     @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component     @param axis_offset: C{axis_offset} to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
2984                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
2985     @type axis_offset: C{int}     @type axis_offset: C{int}
2986     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2987     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
# Line 2939  def trace(arg,axis_offset=0): Line 2989  def trace(arg,axis_offset=0):
2989     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
2990        sh=arg.shape        sh=arg.shape
2991        if len(sh)<2:        if len(sh)<2:
2992          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2993        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
2994          raise ValueError,"trace: axis_offset must be between 0 and %s"%len(sh)-2          raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
2995        s1=1        s1=1
2996        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
2997        s2=1        s2=1
2998        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2999        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
3000          raise ValueError,"trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3001        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3002        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
3003        for i1 in range(s1):        for i1 in range(s1):
# Line 2956  def trace(arg,axis_offset=0): Line 3006  def trace(arg,axis_offset=0):
3006        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3007        return out        return out
3008     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3009        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
3010            raise ValueError,"rank of argument must be greater than 1"
3011          if axis_offset<0 or axis_offset>arg.getRank()-2:
3012            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3013          s=list(arg.getShape())        
3014          if not s[axis_offset] == s[axis_offset+1]:
3015            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3016          return arg._trace(axis_offset)
3017     elif isinstance(arg,float):     elif isinstance(arg,float):
3018        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
3019     elif isinstance(arg,int):     elif isinstance(arg,int):
3020        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
3021     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3022        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
3023     else:     else:
3024        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
3025    
 def escript_trace(arg,axis_offset):  
       "arg is a Data object"  
       if arg.getRank()<2:  
         raise ValueError,"escript_trace: rank of argument must be greater than 1"  
       if axis_offset<0 or axis_offset>arg.getRank()-2:  
         raise ValueError,"escript_trace: axis_offset must be between 0 and %s"%arg.getRank()-2  
       s=list(arg.getShape())          
       if not s[axis_offset] == s[axis_offset+1]:  
         raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)  
       out=arg._matrixtrace(axis_offset)  
       return out  
3026  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
3027     """     """
3028     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 2986  class Trace_Symbol(DependendSymbol): Line 3032  class Trace_Symbol(DependendSymbol):
3032        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3033        @param arg: argument of function        @param arg: argument of function
3034        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3035        @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component        @param axis_offset: C{axis_offset} to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
3036                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3037        @type axis_offset: C{int}        @type axis_offset: C{int}
3038        """        """
3039        if arg.getRank()<2:        if arg.getRank()<2:
3040          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3041        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3042          raise ValueError,"Trace_Symbol: axis_offset must be between 0 and %s"%arg.getRank()-2          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3043        s=list(arg.getShape())                s=list(arg.getShape())        
3044        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3045          raise ValueError,"Trace_Symbol: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3046        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3047    
3048     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3053  class Trace_Symbol(DependendSymbol): Line 3099  class Trace_Symbol(DependendSymbol):
3099    
3100  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3101     """     """
3102     returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.     returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3103    
3104     @param arg: argument     @param arg: argument
3105     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3106     @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.     @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3107                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3108     @type axis_offset: C{int}     @type axis_offset: C{int}
3109     @return: transpose of arg     @return: transpose of arg
3110     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
# Line 3067  def transpose(arg,axis_offset=None): Line 3113  def transpose(arg,axis_offset=None):
3113        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3114        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3115     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3116        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3117        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3118          if axis_offset<0 or axis_offset>r:
3119            raise ValueError,"axis_offset must be between 0 and %s"%r
3120          return arg._transpose(axis_offset)
3121     elif isinstance(arg,float):     elif isinstance(arg,float):
3122        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3123          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3124        return arg        return arg
3125     elif isinstance(arg,int):     elif isinstance(arg,int):
3126        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3127          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3128        return float(arg)        return float(arg)
3129     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3130        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3131        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3132     else:     else:
3133        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3134    
 def escript_transpose(arg,axis_offset):  
       "arg is a Data object"  
       r=arg.getRank()  
       if axis_offset<0 or axis_offset>r:  
         raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r  
       out=arg._transpose(axis_offset)  
       return out  
3135  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3136     """     """
3137     L{Symbol} representing the result of the transpose function     L{Symbol} representing the result of the transpose function
# Line 3100  class Transpose_Symbol(DependendSymbol): Line 3142  class Transpose_Symbol(DependendSymbol):
3142    
3143        @param arg: argument of function        @param arg: argument of function
3144        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3145        @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.        @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3146                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3147        @type axis_offset: C{int}        @type axis_offset: C{int}
3148        """        """
3149        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3150        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3151          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3152        s=arg.getShape()        s=arg.getShape()
3153        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3154    
# Line 3161  class Transpose_Symbol(DependendSymbol): Line 3203  class Transpose_Symbol(DependendSymbol):
3203           return identity(self.getShape())           return identity(self.getShape())
3204        else:        else:
3205           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3206    
3207    def swap_axes(arg,axis0=0,axis1=1):
3208       """
3209       returns the swap of arg by swaping the components axis0 and axis1
3210    
3211       @param arg: argument
3212       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3213       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3214       @type axis0: C{int}
3215       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3216       @type axis1: C{int}
3217       @return: C{arg} with swaped components
3218       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3219       """
3220       if axis0 > axis1:
3221          axis0,axis1=axis1,axis0
3222       if isinstance(arg,numarray.NumArray):
3223          return numarray.swapaxes(arg,axis0,axis1)
3224       elif isinstance(arg,escript.Data):
3225          return arg._swap_axes(axis0,axis1)
3226       elif isinstance(arg,float):
3227          raise TyepError,"float argument is not supported."
3228       elif isinstance(arg,int):
3229          raise TyepError,"int argument is not supported."
3230       elif isinstance(arg,Symbol):
3231          return SwapAxes_Symbol(arg,axis0,axis1)
3232       else:
3233          raise TypeError,"Unknown argument type."
3234    
3235    class SwapAxes_Symbol(DependendSymbol):
3236       """
3237       L{Symbol} representing the result of the swap function
3238       """
3239       def __init__(self,arg,axis0=0,axis1=1):
3240          """
3241          initialization of swap L{Symbol} with argument arg
3242    
3243          @param arg: argument
3244          @type arg: L{Symbol}.
3245          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3246          @type axis0: C{int}
3247          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3248          @type axis1: C{int}
3249          """
3250          if arg.getRank()<2:
3251             raise ValueError,"argument must have at least rank 2."
3252          if axis0<0 or axis0>arg.getRank()-1:
3253             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3254          if axis1<0 or axis1>arg.getRank()-1:
3255             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3256          if axis0 == axis1:
3257             raise ValueError,"axis indices must be different."
3258          if axis0 > axis1:
3259             axis0,axis1=axis1,axis0
3260          s=arg.getShape()
3261          s_out=[]
3262          for i in range(len(s)):
3263             if i == axis0:
3264                s_out.append(s[axis1])
3265             elif i == axis1:
3266                s_out.append(s[axis0])
3267             else:
3268                s_out.append(s[i])
3269          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3270    
3271       def getMyCode(self,argstrs,format="escript"):
3272          """
3273          returns a program code that can be used to evaluate the symbol.
3274    
3275          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3276          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3277          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3278          @type format: C{str}
3279          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3280          @rtype: C{str}
3281          @raise NotImplementedError: if the requested format is not available
3282          """
3283          if format=="escript" or format=="str"  or format=="text":
3284             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3285          else:
3286             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3287    
3288       def substitute(self,argvals):
3289          """
3290          assigns new values to symbols in the definition of the symbol.
3291          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3292    
3293          @param argvals: new values assigned to symbols
3294          @type argvals: C{dict} with keywords of type L{Symbol}.
3295          @return: result of the substitution process. Operations are executed as much as possible.
3296          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3297          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3298          """
3299          if argvals.has_key(self):
3300             arg=argvals[self]
3301             if self.isAppropriateValue(arg):
3302                return arg
3303             else:
3304                raise TypeError,"%s: new value is not appropriate."%str(self)
3305          else:
3306             arg=self.getSubstitutedArguments(argvals)
3307             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3308    
3309       def diff(self,arg):
3310          """
3311          differential of this object
3312    
3313          @param arg: the derivative is calculated with respect to arg
3314          @type arg: L{escript.Symbol}
3315          @return: derivative with respect to C{arg}
3316          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3317          """
3318          if arg==self:
3319             return identity(self.getShape())
3320          else:
3321             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3322    
3323  def symmetric(arg):  def symmetric(arg):
3324      """      """
3325      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
# Line 3173  def symmetric(arg): Line 3332  def symmetric(arg):
3332      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3333        if arg.rank==2:        if arg.rank==2:
3334          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3335             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3336        elif arg.rank==4:        elif arg.rank==4:
3337          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3338             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3339        else:        else:
3340          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3341        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3342      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3343        return escript_symmetric(arg)        if arg.getRank()==2:
3344            if not (arg.getShape()[0]==arg.getShape()[1]):
3345               raise ValueError,"argument must be square."
3346            return arg._symmetric()
3347          elif arg.getRank()==4:
3348            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3349               raise ValueError,"argument must be square."
3350            return arg._symmetric()
3351          else:
3352            raise ValueError,"rank 2 or 4 is required."
3353      elif isinstance(arg,float):      elif isinstance(arg,float):
3354        return arg        return arg
3355      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3189  def symmetric(arg): Line 3357  def symmetric(arg):
3357      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3358        if arg.getRank()==2:        if arg.getRank()==2:
3359          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3360             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3361        elif arg.getRank()==4:        elif arg.getRank()==4:
3362          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3363             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3364        else:        else:
3365          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3366        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3367      else:      else:
3368        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3369    
 def escript_symmetric(arg):  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=arg._symmetric()  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=arg._symmetric()  
       else:  
         raise ValueError,"escript_symmetric: rank 2 or 4 is required."  
       return out  
   
3370  def nonsymmetric(arg):  def nonsymmetric(arg):
3371      """      """
3372      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
# Line 3232  def nonsymmetric(arg): Line 3387  def nonsymmetric(arg):
3387          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3388        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3389      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3390        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3391            if not (arg.getShape()[0]==arg.getShape()[1]):
3392               raise ValueError,"argument must be square."
3393            return arg._nonsymmetric()
3394          elif arg.getRank()==4:
3395            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3396               raise ValueError,"argument must be square."
3397            return arg._nonsymmetric()
3398          else:
3399            raise ValueError,"rank 2 or 4 is required."
3400      elif isinstance(arg,float):      elif isinstance(arg,float):
3401        return arg        return arg
3402      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3250  def nonsymmetric(arg): Line 3414  def nonsymmetric(arg):
3414      else:      else:
3415        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3416    
 def escript_nonsymmetric(arg): # this should be implemented in c++  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=arg._nonsymmetric()  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=arg._nonsymmetric()  
       else:  
         raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."  
       return out  
   
   
3417  def inverse(arg):  def inverse(arg):
3418      """      """
3419      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3420    
3421      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3422      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3423      @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])      @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3424      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3425      @note: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3426      """      """
# Line 3407  class Inverse_Symbol(DependendSymbol): Line 3557  class Inverse_Symbol(DependendSymbol):
3557        if arg==self:        if arg==self:
3558           return identity(self.getShape())           return identity(self.getShape())
3559        else:        else:
3560           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3561    
3562  def eigenvalues(arg):  def eigenvalues(arg):
3563      """      """
# Line 3545  class Add_Symbol(DependendSymbol): Line 3695  class Add_Symbol(DependendSymbol):
3695         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3696         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3697         """         """
3698         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3699         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3700         if not sh0==sh1:         if not sh0==sh1:
3701            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3702         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3620  def mult(arg0,arg1): Line 3770  def mult(arg0,arg1):
3770         """         """
3771         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3772         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3773            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3774         else:         else:
3775            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3776                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3644  class Mult_Symbol(DependendSymbol): Line 3794  class Mult_Symbol(DependendSymbol):
3794         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3795         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3796         """         """
3797         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3798         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3799         if not sh0==sh1:         if not sh0==sh1:
3800            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3801         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3720  def quotient(arg0,arg1): Line 3870  def quotient(arg0,arg1):
3870         """         """
3871         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3872         if testForZero(args[0]):         if testForZero(args[0]):
3873            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3874         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3875            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3876               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3749  class Quotient_Symbol(DependendSymbol): Line 3899  class Quotient_Symbol(DependendSymbol):
3899         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3900         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3901         """         """
3902         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3903         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3904         if not sh0==sh1:         if not sh0==sh1:
3905            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3906         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3826  def power(arg0,arg1): Line 3976  def power(arg0,arg1):
3976         """         """
3977         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3978         if testForZero(args[0]):         if testForZero(args[0]):
3979            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3980         elif testForZero(args[1]):         elif testForZero(args[1]):
3981            return numarray.ones(pokeShape(args[1]),numarray.Float64)            return numarray.ones(getShape(args[1]),numarray.Float64)
3982         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3983            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3984         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
# Line 3851  class Power_Symbol(DependendSymbol): Line 4001  class Power_Symbol(DependendSymbol):
4001         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4002         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4003         """         """
4004         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4005         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4006         if not sh0==sh1:         if not sh0==sh1:
4007            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4008         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3951  def minimum(*args): Line 4101  def minimum(*args):
4101            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4102      return out      return out
4103    
4104  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4105      """      """
4106      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4107    
4108      @param arg: argument      @param arg: argument
4109      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4110      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4111      @type minval: C{float}      @type minval: C{float} or C{None}
4112      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4113      @type maxval: C{float}      @type maxval: C{float} or C{None}
4114      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4115               less then maxval are unchanged.               less then maxval are unchanged.
4116      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4117      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4118      """      """
4119      if minval>maxval:      if not minval==None and not maxval==None:
4120         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4121      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4122        if minval == None:
4123            tmp=arg
4124        else:
4125            tmp=maximum(minval,arg)
4126        if maxval == None:
4127            return tmp
4128        else:
4129            return minimum(tmp,maxval)
4130    
4131        
4132  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3989  def inner(arg0,arg1): Line 4147  def inner(arg0,arg1):
4147      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input
4148      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4149      """      """
4150      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4151      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4152      if not sh0==sh1:      if not sh0==sh1:
4153          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4154      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4155    
4156    def outer(arg0,arg1):
4157        """
4158        the outer product of the two argument:
4159    
4160        out[t,s]=arg0[t]*arg1[s]
4161    
4162        where
4163    
4164            - s runs through arg0.Shape
4165            - t runs through arg1.Shape
4166    
4167        @param arg0: first argument
4168        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4169        @param arg1: second argument
4170        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4171        @return: the outer product of arg0 and arg1 at each data point
4172        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4173        """
4174        return generalTensorProduct(arg0,arg1,axis_offset=0)
4175    
4176  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4177      """      """
4178        see L{matrix_mult}
4179        """
4180        return matrix_mult(arg0,arg1)
4181    
4182    def matrix_mult(arg0,arg1):
4183        """
4184      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4185    
4186      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
# Line 4005  def matrixmult(arg0,arg1): Line 4189  def matrixmult(arg0,arg1):
4189    
4190      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4191    
4192      The second dimension of arg0 and the length of arg1 must match.      The second dimension of arg0 and the first dimension of arg1 must match.
4193    
4194      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4195      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4015  def matrixmult(arg0,arg1): Line 4199  def matrixmult(arg0,arg1):
4199      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4200      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4201      """      """
4202      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4203      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4204      if not len(sh0)==2 :      if not len(sh0)==2 :
4205          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4206      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4207          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4208      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4209    
4210  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4211      """      """
4212      the outer product of the two argument:      use L{tensor_mult}
   
     out[t,s]=arg0[t]*arg1[s]  
   
     where  
   
         - s runs through arg0.Shape  
         - t runs through arg1.Shape  
   
     @param arg0: first argument  
     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}  
     @param arg1: second argument  
     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}  
     @return: the outer product of arg0 and arg1 at each data point  
     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input  
4213      """      """
4214      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
4215    
4216    def tensor_mult(arg0,arg1):
 def tensormult(arg0,arg1):  
4217      """      """
4218      the tensor product of the two argument:      the tensor product of the two argument:
4219            
# Line 4069  def tensormult(arg0,arg1): Line 4238  def tensormult(arg0,arg1):
4238    
4239      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]
4240    
4241      In the first case the the second dimension of arg0 and the length of arg1 must match and        In the first case the the second dimension of arg0 and the last dimension of arg1 must match and  
4242      in the second case the two last dimensions of arg0 must match the shape of arg1.      in the second case the two last dimensions of arg0 must match the two first dimensions of arg1.
4243    
4244      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4245      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4079  def tensormult(arg0,arg1): Line 4248  def tensormult(arg0,arg1):
4248      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4249      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4250      """      """
4251      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4252      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4253      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4254         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4255      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4256         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4257      else:      else:
4258          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4259    
4260  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4261      """      """
# Line 4100  def generalTensorProduct(arg0,arg1,axis_ Line 4269  def generalTensorProduct(arg0,arg1,axis_
4269          - r runs trough arg0.Shape[:axis_offset]          - r runs trough arg0.Shape[:axis_offset]
4270          - t runs through arg1.Shape[axis_offset:]          - t runs through arg1.Shape[axis_offset:]
4271    
     In the first case the the second dimension of arg0 and the length of arg1 must match and    
     in the second case the two last dimensions of arg0 must match the shape of arg1.  
   
4272      @param arg0: first argument      @param arg0: first argument
4273      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4274      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4275      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4276      @return: the general tensor product of arg0 and arg1 at each data point.      @return: the general tensor product of arg0 and arg1 at each data point.
4277      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
# Line 4118  def generalTensorProduct(arg0,arg1,axis_ Line 4284  def generalTensorProduct(arg0,arg1,axis_
4284             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4285         else:         else:
4286             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4287                 raise ValueError,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)                 raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4288             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4289             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4290             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4144  def generalTensorProduct(arg0,arg1,axis_ Line 4310  def generalTensorProduct(arg0,arg1,axis_
4310                                    
4311  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4312     """     """
4313     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4314     """     """
4315     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4316         """         """
4317         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4318    
4319         @param arg0: numerator         @param arg0: first argument
4320         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4321         @param arg1: denominator         @param arg1: second argument
4322         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4323         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4324         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4325         """         """
4326         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4327         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4328         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4329         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4330         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 4205  class GeneralTensorProduct_Symbol(Depend Line 4371  class GeneralTensorProduct_Symbol(Depend
4371           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4372           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4373    
4374  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4375      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4376      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4377      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4378      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4379      shape10=arg1.getShape()[:axis_offset]      """
4380      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4381      if not shape01==shape10:  
4382          raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4383    
4384      # whatr function space should be used? (this here is not good!)      or
4385      fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
4386      # create return value:      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4387      out=escript.Data(0.,tuple(shape0+shape1),fs)  
4388      #      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4389      s0=[[]]  
4390      for k in shape0:      The first dimension of arg0 and arg1 must match.
4391            s=[]  
4392            for j in s0:      @param arg0: first argument of rank 2
4393                  for i in range(k): s.append(j+[slice(i,i)])      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4394            s0=s      @param arg1: second argument of at least rank 1
4395      s1=[[]]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4396      for k in shape1:      @return: the product of the transposed of arg0 and arg1 at each data point
4397            s=[]      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4398            for j in s1:      @raise ValueError: if the shapes of the arguments are not appropriate
4399                  for i in range(k): s.append(j+[slice(i,i)])      """
4400            s1=s      sh0=getShape(arg0)
4401      s01=[[]]      sh1=getShape(arg1)
4402      for k in shape01:      if not len(sh0)==2 :
4403            s=[]          raise ValueError,"first argument must have rank 2"
4404            for j in s01:      if not len(sh1)==2 and not len(sh1)==1:
4405                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"second argument must have rank 1 or 2"
4406            s01=s      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4407    
4408      for i0 in s0:  def transposed_tensor_mult(arg0,arg1):
4409         for i1 in s1:      """
4410           s=escript.Scalar(0.,fs)      the tensor product of the transposed of the first and the second argument
4411           for i01 in s01:      
4412              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      for arg0 of rank 2 this is
4413           out.__setitem__(tuple(i0+i1),s)  
4414      return out      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4415    
4416        or
4417    
4418        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4419    
4420      
4421        and for arg0 of rank 4 this is
4422    
4423        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4424                  
4425        or
4426    
4427        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4428    
4429        or
4430    
4431        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4432    
4433        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4434        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4435    
4436        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4437    
4438        @param arg0: first argument of rank 2 or 4
4439        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4440        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4441        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4442        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4443        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4444        """
4445        sh0=getShape(arg0)
4446        sh1=getShape(arg1)
4447        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4448           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4449        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4450           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4451        else:
4452            raise ValueError,"first argument must have rank 2 or 4"
4453    
4454    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4455        """
4456        generalized tensor product of transposed of arg0 and arg1:
4457    
4458        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4459    
4460        where
4461    
4462            - s runs through arg0.Shape[axis_offset:]
4463            - r runs trough arg0.Shape[:axis_offset]
4464            - t runs through arg1.Shape[axis_offset:]
4465    
4466        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4467        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4468    
4469        @param arg0: first argument
4470        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4471        @param arg1: second argument
4472        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4473        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4474        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4475        """
4476        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4477        arg0,arg1=matchType(arg0,arg1)
4478        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4479        if isinstance(arg0,numarray.NumArray):
4480           if isinstance(arg1,Symbol):
4481               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4482           else:
4483               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4484                   raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4485               arg0_c=arg0.copy()
4486               arg1_c=arg1.copy()
4487               sh0,sh1=arg0.shape,arg1.shape
4488               d0,d1,d01=1,1,1
4489               for i in sh0[axis_offset:]: d0*=i
4490               for i in sh1[axis_offset:]: d1*=i
4491               for i in sh0[:axis_offset]: d01*=i
4492               arg0_c.resize((d01,d0))
4493               arg1_c.resize((d01,d1))
4494               out=numarray.zeros((d0,d1),numarray.Float64)
4495               for i0 in range(d0):
4496                        for i1 in range(d1):
4497                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4498               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4499               return out
4500        elif isinstance(arg0,escript.Data):
4501           if isinstance(arg1,Symbol):
4502               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4503           else:
4504               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4505        else:      
4506           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4507                    
4508    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4509       """
4510       Symbol representing the general tensor product of the transposed of arg0 and arg1
4511       """
4512       def __init__(self,arg0,arg1,axis_offset=0):
4513           """
4514           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4515    
4516           @param arg0: first argument
4517           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4518           @param arg1: second argument
4519           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4520           @raise ValueError: inconsistent dimensions of arguments.
4521           @note: if both arguments have a spatial dimension, they must equal.
4522           """
4523           sh_arg0=getShape(arg0)
4524           sh_arg1=getShape(arg1)
4525           sh01=sh_arg0[:axis_offset]
4526           sh10=sh_arg1[:axis_offset]
4527           sh0=sh_arg0[axis_offset:]
4528           sh1=sh_arg1[axis_offset:]
4529           if not sh01==sh10:
4530               raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4531           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4532    
4533       def getMyCode(self,argstrs,format="escript"):
4534          """
4535          returns a program code that can be used to evaluate the symbol.
4536    
4537          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4538          @type argstrs: C{list} of length 2 of C{str}.
4539          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4540          @type format: C{str}
4541          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4542          @rtype: C{str}
4543          @raise NotImplementedError: if the requested format is not available
4544          """
4545          if format=="escript" or format=="str" or format=="text":
4546             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4547          else:
4548             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4549    
4550       def substitute(self,argvals):
4551          """
4552          assigns new values to symbols in the definition of the symbol.
4553          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4554    
4555          @param argvals: new values assigned to symbols
4556          @type argvals: C{dict} with keywords of type L{Symbol}.
4557          @return: result of the substitution process. Operations are executed as much as possible.
4558          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4559          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4560          """
4561          if argvals.has_key(self):
4562             arg=argvals[self]
4563             if self.isAppropriateValue(arg):
4564                return arg
4565             else:
4566                raise TypeError,"%s: new value is not appropriate."%str(self)
4567          else:
4568             args=self.getSubstitutedArguments(argvals)
4569             return generalTransposedTensorProduct(args[0],args[1],args[2])
4570    
4571    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4572        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4573        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4574    
4575    def matrix_transposed_mult(arg0,arg1):
4576        """
4577        matrix-transposed(matrix) product of the two argument:
4578    
4579        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4580    
4581        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4582    
4583        The last dimensions of arg0 and arg1 must match.
4584    
4585        @param arg0: first argument of rank 2
4586        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4587        @param arg1: second argument of rank 2
4588        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4589        @return: the product of arg0 and the transposed of arg1 at each data point
4590        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4591        @raise ValueError: if the shapes of the arguments are not appropriate
4592        """
4593        sh0=getShape(arg0)
4594        sh1=getShape(arg1)
4595        if not len(sh0)==2 :
4596            raise ValueError,"first argument must have rank 2"
4597        if not len(sh1)==2 and not len(sh1)==1:
4598            raise ValueError,"second argument must have rank 1 or 2"
4599        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4600    
4601    def tensor_transposed_mult(arg0,arg1):
4602        """
4603        the tensor product of the first and the transpose of the second argument
4604        
4605        for arg0 of rank 2 this is
4606    
4607        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4608    
4609        and for arg0 of rank 4 this is
4610    
4611        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4612                  
4613        or
4614    
4615        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4616    
4617        In the first case the the second dimension of arg0 and arg1 must match and  
4618        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4619    
4620        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4621    
4622        @param arg0: first argument of rank 2 or 4
4623        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4624        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4625        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4626        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4627        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4628        """
4629        sh0=getShape(arg0)
4630        sh1=getShape(arg1)
4631        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4632           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4633        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4634           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4635        else:
4636            raise ValueError,"first argument must have rank 2 or 4"
4637    
4638    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4639        """
4640        generalized tensor product of transposed of arg0 and arg1:
4641    
4642        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4643    
4644        where
4645    
4646            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4647            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4648            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4649    
4650        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4651        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4652    
4653        @param arg0: first argument
4654        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4655        @param arg1: second argument
4656        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4657        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4658        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4659        """
4660        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4661        arg0,arg1=matchType(arg0,arg1)
4662        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4663        if isinstance(arg0,numarray.NumArray):
4664           if isinstance(arg1,Symbol):
4665               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4666           else:
4667               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4668                   raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4669               arg0_c=arg0.copy()
4670               arg1_c=arg1.copy()
4671               sh0,sh1=arg0.shape,arg1.shape
4672               d0,d1,d01=1,1,1
4673               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4674               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4675               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4676               arg0_c.resize((d0,d01))
4677               arg1_c.resize((d1,d01))
4678               out=numarray.zeros((d0,d1),numarray.Float64)
4679               for i0 in range(d0):
4680                        for i1 in range(d1):
4681                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4682               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4683               return out
4684        elif isinstance(arg0,escript.Data):
4685           if isinstance(arg1,Symbol):
4686               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4687           else:
4688               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4689        else:      
4690           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4691                    
4692    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4693       """
4694       Symbol representing the general tensor product of arg0 and the transpose of arg1
4695       """
4696       def __init__(self,arg0,arg1,axis_offset=0):
4697           """
4698           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4699    
4700           @param arg0: first argument
4701           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4702           @param arg1: second argument
4703           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4704           @raise ValueError: inconsistent dimensions of arguments.
4705           @note: if both arguments have a spatial dimension, they must equal.
4706           """
4707           sh_arg0=getShape(arg0)
4708           sh_arg1=getShape(arg1)
4709           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4710           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4711           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4712           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4713           if not sh01==sh10:
4714               raise ValueError,"dimensions of last %s components in left argument don't match the last %s components in the right argument."%(axis_offset,axis_offset)
4715           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4716    
4717       def getMyCode(self,argstrs,format="escript"):
4718          """
4719          returns a program code that can be used to evaluate the symbol.
4720    
4721          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4722          @type argstrs: C{list} of length 2 of C{str}.
4723          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4724          @type format: C{str}
4725          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4726          @rtype: C{str}
4727          @raise NotImplementedError: if the requested format is not available
4728          """
4729          if format=="escript" or format=="str" or format=="text":
4730             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4731          else:
4732             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4733    
4734       def substitute(self,argvals):
4735          """
4736          assigns new values to symbols in the definition of the symbol.
4737          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4738    
4739          @param argvals: new values assigned to symbols
4740          @type argvals: C{dict} with keywords of type L{Symbol}.
4741          @return: result of the substitution process. Operations are executed as much as possible.
4742          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4743          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4744          """
4745          if argvals.has_key(self):
4746             arg=argvals[self]
4747             if self.isAppropriateValue(arg):
4748                return arg
4749             else:
4750                raise TypeError,"%s: new value is not appropriate."%str(self)
4751          else:
4752             args=self.getSubstitutedArguments(argvals)
4753             return generalTensorTransposedProduct(args[0],args[1],args[2])
4754    
4755    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4756        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4757        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4758    
4759  #=========================================================  #=========================================================
4760  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency

Legend:
Removed from v.775  
changed lines
  Added in v.1042

  ViewVC Help
Powered by ViewVC 1.1.26