/[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 1044 by gross, Mon Mar 19 07:29:31 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: L{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 insertTagNames(domain,**kwargs):
47        """
48        inserts tag names into the domain
49    
50        @param domain: a domain object
51        @type domain: C{escript.Domain}
52        @keyword <tag name>: tag key assigned to <tag name>
53        @type <tag name>: C{int}
54        """
55        for  k in kwargs:
56             domain.setTagMap(k,kwargs[k])
57    
58    def insertTaggedValues(target,**kwargs):
59        """
60        inserts tagged values into the tagged using tag names
61    
62        @param target: data to be filled by tagged values
63        @type target: L{escript.Data}
64        @keyword <tag name>: value to be used for <tag name>
65        @type <tag name>: C{float} or {numarray.NumArray}
66        @return: C{target}
67        @rtype: L{escript.Data}
68        """
69        for k in kwargs:
70            target.setTaggedValue(k,kwargs[k])
71        return target
72    
73        
74  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
75      """      """
76      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 279  def inf(arg):
279  #=========================================================================  #=========================================================================
280  #   some little helpers  #   some little helpers
281  #=========================================================================  #=========================================================================
282  def pokeShape(arg):  def getRank(arg):
283        """
284        identifies the rank of its argument
285    
286        @param arg: a given object
287        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
288        @return: the rank of the argument
289        @rtype: C{int}
290        @raise TypeError: if type of arg cannot be processed
291        """
292    
293        if isinstance(arg,numarray.NumArray):
294            return arg.rank
295        elif isinstance(arg,escript.Data):
296            return arg.getRank()
297        elif isinstance(arg,float):
298            return 0
299        elif isinstance(arg,int):
300            return 0
301        elif isinstance(arg,Symbol):
302            return arg.getRank()
303        else:
304          raise TypeError,"getShape: cannot identify shape"
305    def getShape(arg):
306      """      """
307      identifies the shape of its argument      identifies the shape of its argument
308    
# Line 260  def pokeShape(arg): Line 324  def pokeShape(arg):
324      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
325          return arg.getShape()          return arg.getShape()
326      else:      else:
327        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
328    
329  def pokeDim(arg):  def pokeDim(arg):
330      """      """
# Line 283  def commonShape(arg0,arg1): Line 347  def commonShape(arg0,arg1):
347      """      """
348      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.
349    
350      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
351      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
352      @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.
353      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
354      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
355      """      """
356      sh0=pokeShape(arg0)      sh0=getShape(arg0)
357      sh1=pokeShape(arg1)      sh1=getShape(arg1)
358      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
359         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
360               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 504  def matchShape(arg0,arg1):
504      @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}
505      @param arg1: a given object      @param arg1: a given object
506      @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}
507      @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.
508      @rtype: C{tuple}      @rtype: C{tuple}
509      """      """
510      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
511      sh0=pokeShape(arg0)      sh0=getShape(arg0)
512      sh1=pokeShape(arg1)      sh1=getShape(arg1)
513      if len(sh0)<len(sh):      if len(sh0)<len(sh):
514         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
515      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
# Line 573  class Symbol(object): Line 637  class Symbol(object):
637            if isinstance(a,Symbol):            if isinstance(a,Symbol):
638               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
639            else:            else:
640                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
641                if len(s)>0:                if len(s)>0:
642                   out.append(numarray.zeros(s),numarray.Float64)                   out.append(numarray.zeros(s),numarray.Float64)
643                else:                else:
# Line 1310  def whereNonZero(arg,tol=0.): Line 1374  def whereNonZero(arg,tol=0.):
1374     else:     else:
1375        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1376    
1377    def erf(arg):
1378       """
1379       returns erf of argument arg
1380    
1381       @param arg: argument
1382       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1383       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1384       @raises TypeError: if the type of the argument is not expected.
1385       """
1386       if isinstance(arg,escript.Data):
1387          return arg._erf()
1388       else:
1389          raise TypeError,"erf: Unknown argument type."
1390    
1391  def sin(arg):  def sin(arg):
1392     """     """
1393     returns sine of argument arg     returns sine of argument arg
# Line 2930  def trace(arg,axis_offset=0): Line 3008  def trace(arg,axis_offset=0):
3008    
3009     @param arg: argument     @param arg: argument
3010     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3011     @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
3012                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3013     @type axis_offset: C{int}     @type axis_offset: C{int}
3014     @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.
3015     @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 3017  def trace(arg,axis_offset=0):
3017     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
3018        sh=arg.shape        sh=arg.shape
3019        if len(sh)<2:        if len(sh)<2:
3020          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3021        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
3022          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
3023        s1=1        s1=1
3024        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
3025        s2=1        s2=1
3026        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3027        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
3028          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)
3029        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))
3030        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
3031        for i1 in range(s1):        for i1 in range(s1):
# Line 2956  def trace(arg,axis_offset=0): Line 3034  def trace(arg,axis_offset=0):
3034        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
3035        return out        return out
3036     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3037        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
3038            raise ValueError,"rank of argument must be greater than 1"
3039          if axis_offset<0 or axis_offset>arg.getRank()-2:
3040            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3041          s=list(arg.getShape())        
3042          if not s[axis_offset] == s[axis_offset+1]:
3043            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3044          return arg._trace(axis_offset)
3045     elif isinstance(arg,float):     elif isinstance(arg,float):
3046        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
3047     elif isinstance(arg,int):     elif isinstance(arg,int):
3048        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
3049     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3050        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
3051     else:     else:
3052        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
3053    
 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  
3054  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
3055     """     """
3056     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 3060  class Trace_Symbol(DependendSymbol):
3060        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3061        @param arg: argument of function        @param arg: argument of function
3062        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3063        @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
3064                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3065        @type axis_offset: C{int}        @type axis_offset: C{int}
3066        """        """
3067        if arg.getRank()<2:        if arg.getRank()<2:
3068          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3069        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3070          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
3071        s=list(arg.getShape())                s=list(arg.getShape())        
3072        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3073          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)
3074        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())
3075    
3076     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3053  class Trace_Symbol(DependendSymbol): Line 3127  class Trace_Symbol(DependendSymbol):
3127    
3128  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3129     """     """
3130     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.
3131    
3132     @param arg: argument     @param arg: argument
3133     @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}
3134     @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.
3135                         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.
3136     @type axis_offset: C{int}     @type axis_offset: C{int}
3137     @return: transpose of arg     @return: transpose of arg
3138     @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 3141  def transpose(arg,axis_offset=None):
3141        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3142        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))
3143     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3144        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3145        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3146          if axis_offset<0 or axis_offset>r:
3147            raise ValueError,"axis_offset must be between 0 and %s"%r
3148          return arg._transpose(axis_offset)
3149     elif isinstance(arg,float):     elif isinstance(arg,float):
3150        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3151          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3152        return arg        return arg
3153     elif isinstance(arg,int):     elif isinstance(arg,int):
3154        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3155          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3156        return float(arg)        return float(arg)
3157     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3158        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3159        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3160     else:     else:
3161        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3162    
 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  
3163  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3164     """     """
3165     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 3170  class Transpose_Symbol(DependendSymbol):
3170    
3171        @param arg: argument of function        @param arg: argument of function
3172        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3173        @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.
3174                         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.
3175        @type axis_offset: C{int}        @type axis_offset: C{int}
3176        """        """
3177        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3178        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3179          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()
3180        s=arg.getShape()        s=arg.getShape()
3181        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())
3182    
# Line 3161  class Transpose_Symbol(DependendSymbol): Line 3231  class Transpose_Symbol(DependendSymbol):
3231           return identity(self.getShape())           return identity(self.getShape())
3232        else:        else:
3233           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3234    
3235    def swap_axes(arg,axis0=0,axis1=1):
3236       """
3237       returns the swap of arg by swaping the components axis0 and axis1
3238    
3239       @param arg: argument
3240       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3241       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3242       @type axis0: C{int}
3243       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3244       @type axis1: C{int}
3245       @return: C{arg} with swaped components
3246       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3247       """
3248       if axis0 > axis1:
3249          axis0,axis1=axis1,axis0
3250       if isinstance(arg,numarray.NumArray):
3251          return numarray.swapaxes(arg,axis0,axis1)
3252       elif isinstance(arg,escript.Data):
3253          return arg._swap_axes(axis0,axis1)
3254       elif isinstance(arg,float):
3255          raise TyepError,"float argument is not supported."
3256       elif isinstance(arg,int):
3257          raise TyepError,"int argument is not supported."
3258       elif isinstance(arg,Symbol):
3259          return SwapAxes_Symbol(arg,axis0,axis1)
3260       else:
3261          raise TypeError,"Unknown argument type."
3262    
3263    class SwapAxes_Symbol(DependendSymbol):
3264       """
3265       L{Symbol} representing the result of the swap function
3266       """
3267       def __init__(self,arg,axis0=0,axis1=1):
3268          """
3269          initialization of swap L{Symbol} with argument arg
3270    
3271          @param arg: argument
3272          @type arg: L{Symbol}.
3273          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3274          @type axis0: C{int}
3275          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3276          @type axis1: C{int}
3277          """
3278          if arg.getRank()<2:
3279             raise ValueError,"argument must have at least rank 2."
3280          if axis0<0 or axis0>arg.getRank()-1:
3281             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3282          if axis1<0 or axis1>arg.getRank()-1:
3283             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3284          if axis0 == axis1:
3285             raise ValueError,"axis indices must be different."
3286          if axis0 > axis1:
3287             axis0,axis1=axis1,axis0
3288          s=arg.getShape()
3289          s_out=[]
3290          for i in range(len(s)):
3291             if i == axis0:
3292                s_out.append(s[axis1])
3293             elif i == axis1:
3294                s_out.append(s[axis0])
3295             else:
3296                s_out.append(s[i])
3297          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3298    
3299       def getMyCode(self,argstrs,format="escript"):
3300          """
3301          returns a program code that can be used to evaluate the symbol.
3302    
3303          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3304          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3305          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3306          @type format: C{str}
3307          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3308          @rtype: C{str}
3309          @raise NotImplementedError: if the requested format is not available
3310          """
3311          if format=="escript" or format=="str"  or format=="text":
3312             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3313          else:
3314             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3315    
3316       def substitute(self,argvals):
3317          """
3318          assigns new values to symbols in the definition of the symbol.
3319          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3320    
3321          @param argvals: new values assigned to symbols
3322          @type argvals: C{dict} with keywords of type L{Symbol}.
3323          @return: result of the substitution process. Operations are executed as much as possible.
3324          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3325          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3326          """
3327          if argvals.has_key(self):
3328             arg=argvals[self]
3329             if self.isAppropriateValue(arg):
3330                return arg
3331             else:
3332                raise TypeError,"%s: new value is not appropriate."%str(self)
3333          else:
3334             arg=self.getSubstitutedArguments(argvals)
3335             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3336    
3337       def diff(self,arg):
3338          """
3339          differential of this object
3340    
3341          @param arg: the derivative is calculated with respect to arg
3342          @type arg: L{escript.Symbol}
3343          @return: derivative with respect to C{arg}
3344          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3345          """
3346          if arg==self:
3347             return identity(self.getShape())
3348          else:
3349             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3350    
3351  def symmetric(arg):  def symmetric(arg):
3352      """      """
3353      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 3360  def symmetric(arg):
3360      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3361        if arg.rank==2:        if arg.rank==2:
3362          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3363             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3364        elif arg.rank==4:        elif arg.rank==4:
3365          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]):
3366             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3367        else:        else:
3368          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3369        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3370      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3371        return escript_symmetric(arg)        if arg.getRank()==2:
3372            if not (arg.getShape()[0]==arg.getShape()[1]):
3373               raise ValueError,"argument must be square."
3374            return arg._symmetric()
3375          elif arg.getRank()==4:
3376            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3377               raise ValueError,"argument must be square."
3378            return arg._symmetric()
3379          else:
3380            raise ValueError,"rank 2 or 4 is required."
3381      elif isinstance(arg,float):      elif isinstance(arg,float):
3382        return arg        return arg
3383      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3189  def symmetric(arg): Line 3385  def symmetric(arg):
3385      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3386        if arg.getRank()==2:        if arg.getRank()==2:
3387          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3388             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3389        elif arg.getRank()==4:        elif arg.getRank()==4:
3390          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]):
3391             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3392        else:        else:
3393          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3394        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3395      else:      else:
3396        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3397    
 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  
   
3398  def nonsymmetric(arg):  def nonsymmetric(arg):
3399      """      """
3400      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 3415  def nonsymmetric(arg):
3415          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3416        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3417      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3418        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3419            if not (arg.getShape()[0]==arg.getShape()[1]):
3420               raise ValueError,"argument must be square."
3421            return arg._nonsymmetric()
3422          elif arg.getRank()==4:
3423            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3424               raise ValueError,"argument must be square."
3425            return arg._nonsymmetric()
3426          else:
3427            raise ValueError,"rank 2 or 4 is required."
3428      elif isinstance(arg,float):      elif isinstance(arg,float):
3429        return arg        return arg
3430      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3250  def nonsymmetric(arg): Line 3442  def nonsymmetric(arg):
3442      else:      else:
3443        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3444    
 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  
   
   
3445  def inverse(arg):  def inverse(arg):
3446      """      """
3447      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3448    
3449      @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.
3450      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3451      @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])
3452      @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
3453      @note: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3454      """      """
# Line 3407  class Inverse_Symbol(DependendSymbol): Line 3585  class Inverse_Symbol(DependendSymbol):
3585        if arg==self:        if arg==self:
3586           return identity(self.getShape())           return identity(self.getShape())
3587        else:        else:
3588           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3589    
3590  def eigenvalues(arg):  def eigenvalues(arg):
3591      """      """
# Line 3545  class Add_Symbol(DependendSymbol): Line 3723  class Add_Symbol(DependendSymbol):
3723         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3724         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3725         """         """
3726         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3727         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3728         if not sh0==sh1:         if not sh0==sh1:
3729            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3730         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 3798  def mult(arg0,arg1):
3798         """         """
3799         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3800         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3801            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3802         else:         else:
3803            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3804                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3644  class Mult_Symbol(DependendSymbol): Line 3822  class Mult_Symbol(DependendSymbol):
3822         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3823         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3824         """         """
3825         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3826         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3827         if not sh0==sh1:         if not sh0==sh1:
3828            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3829         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 3898  def quotient(arg0,arg1):
3898         """         """
3899         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3900         if testForZero(args[0]):         if testForZero(args[0]):
3901            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3902         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3903            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3904               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3749  class Quotient_Symbol(DependendSymbol): Line 3927  class Quotient_Symbol(DependendSymbol):
3927         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3928         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3929         """         """
3930         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3931         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3932         if not sh0==sh1:         if not sh0==sh1:
3933            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3934         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 4004  def power(arg0,arg1):
4004         """         """
4005         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
4006         if testForZero(args[0]):         if testForZero(args[0]):
4007            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
4008         elif testForZero(args[1]):         elif testForZero(args[1]):
4009            return numarray.ones(pokeShape(args[1]),numarray.Float64)            return numarray.ones(getShape(args[1]),numarray.Float64)
4010         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
4011            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
4012         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 4029  class Power_Symbol(DependendSymbol):
4029         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4030         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4031         """         """
4032         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4033         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4034         if not sh0==sh1:         if not sh0==sh1:
4035            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4036         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3951  def minimum(*args): Line 4129  def minimum(*args):
4129            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4130      return out      return out
4131    
4132  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4133      """      """
4134      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4135    
4136      @param arg: argument      @param arg: argument
4137      @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}
4138      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4139      @type minval: C{float}      @type minval: C{float} or C{None}
4140      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4141      @type maxval: C{float}      @type maxval: C{float} or C{None}
4142      @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
4143               less then maxval are unchanged.               less then maxval are unchanged.
4144      @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
4145      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4146      """      """
4147      if minval>maxval:      if not minval==None and not maxval==None:
4148         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4149      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4150        if minval == None:
4151            tmp=arg
4152        else:
4153            tmp=maximum(minval,arg)
4154        if maxval == None:
4155            return tmp
4156        else:
4157            return minimum(tmp,maxval)
4158    
4159        
4160  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3989  def inner(arg0,arg1): Line 4175  def inner(arg0,arg1):
4175      @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
4176      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4177      """      """
4178      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4179      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4180      if not sh0==sh1:      if not sh0==sh1:
4181          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4182      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4183    
4184    def outer(arg0,arg1):
4185        """
4186        the outer product of the two argument:
4187    
4188        out[t,s]=arg0[t]*arg1[s]
4189    
4190        where
4191    
4192            - s runs through arg0.Shape
4193            - t runs through arg1.Shape
4194    
4195        @param arg0: first argument
4196        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4197        @param arg1: second argument
4198        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4199        @return: the outer product of arg0 and arg1 at each data point
4200        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4201        """
4202        return generalTensorProduct(arg0,arg1,axis_offset=0)
4203    
4204  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4205      """      """
4206        see L{matrix_mult}
4207        """
4208        return matrix_mult(arg0,arg1)
4209    
4210    def matrix_mult(arg0,arg1):
4211        """
4212      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4213    
4214      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 4217  def matrixmult(arg0,arg1):
4217    
4218      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4219    
4220      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.
4221    
4222      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4223      @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 4227  def matrixmult(arg0,arg1):
4227      @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
4228      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4229      """      """
4230      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4231      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4232      if not len(sh0)==2 :      if not len(sh0)==2 :
4233          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4234      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
4235          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4236      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4237    
4238  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4239      """      """
4240      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  
4241      """      """
4242      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
   
4243    
4244  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4245      """      """
4246      the tensor product of the two argument:      the tensor product of the two argument:
4247            
# Line 4069  def tensormult(arg0,arg1): Line 4266  def tensormult(arg0,arg1):
4266    
4267      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]
4268    
4269      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  
4270      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.
4271    
4272      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4273      @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 4276  def tensormult(arg0,arg1):
4276      @return: the tensor product of arg0 and arg1 at each data point      @return: the 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
4278      """      """
4279      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4280      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4281      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4282         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4283      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):
4284         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4285      else:      else:
4286          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4287    
4288  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4289      """      """
# Line 4100  def generalTensorProduct(arg0,arg1,axis_ Line 4297  def generalTensorProduct(arg0,arg1,axis_
4297          - r runs trough arg0.Shape[:axis_offset]          - r runs trough arg0.Shape[:axis_offset]
4298          - t runs through arg1.Shape[axis_offset:]          - t runs through arg1.Shape[axis_offset:]
4299    
     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.  
   
4300      @param arg0: first argument      @param arg0: first argument
4301      @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}
4302      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4303      @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}
4304      @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.
4305      @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 4312  def generalTensorProduct(arg0,arg1,axis_
4312             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4313         else:         else:
4314             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4315                 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)
4316             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4317             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4318             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4144  def generalTensorProduct(arg0,arg1,axis_ Line 4338  def generalTensorProduct(arg0,arg1,axis_
4338                                    
4339  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4340     """     """
4341     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4342     """     """
4343     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4344         """         """
4345         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4346    
4347         @param arg0: numerator         @param arg0: first argument
4348         @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}.
4349         @param arg1: denominator         @param arg1: second argument
4350         @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}.
4351         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4352         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4353         """         """
4354         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4355         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4356         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4357         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4358         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 4205  class GeneralTensorProduct_Symbol(Depend Line 4399  class GeneralTensorProduct_Symbol(Depend
4399           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4400           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4401    
4402  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4403      "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!!!"
4404      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4405      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4406      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4407      shape10=arg1.getShape()[:axis_offset]      """
4408      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4409      if not shape01==shape10:  
4410          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]
4411    
4412      # whatr function space should be used? (this here is not good!)      or
4413      fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
4414      # create return value:      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4415      out=escript.Data(0.,tuple(shape0+shape1),fs)  
4416      #      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4417      s0=[[]]  
4418      for k in shape0:      The first dimension of arg0 and arg1 must match.
4419            s=[]  
4420            for j in s0:      @param arg0: first argument of rank 2
4421                  for i in range(k): s.append(j+[slice(i,i)])      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4422            s0=s      @param arg1: second argument of at least rank 1
4423      s1=[[]]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4424      for k in shape1:      @return: the product of the transposed of arg0 and arg1 at each data point
4425            s=[]      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4426            for j in s1:      @raise ValueError: if the shapes of the arguments are not appropriate
4427                  for i in range(k): s.append(j+[slice(i,i)])      """
4428            s1=s      sh0=getShape(arg0)
4429      s01=[[]]      sh1=getShape(arg1)
4430      for k in shape01:      if not len(sh0)==2 :
4431            s=[]          raise ValueError,"first argument must have rank 2"
4432            for j in s01:      if not len(sh1)==2 and not len(sh1)==1:
4433                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"second argument must have rank 1 or 2"
4434            s01=s      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4435    
4436      for i0 in s0:  def transposed_tensor_mult(arg0,arg1):
4437         for i1 in s1:      """
4438           s=escript.Scalar(0.,fs)      the tensor product of the transposed of the first and the second argument
4439           for i01 in s01:      
4440              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      for arg0 of rank 2 this is
4441           out.__setitem__(tuple(i0+i1),s)  
4442      return out      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4443    
4444        or
4445    
4446        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4447    
4448      
4449        and for arg0 of rank 4 this is
4450    
4451        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4452                  
4453        or
4454    
4455        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4456    
4457        or
4458    
4459        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4460    
4461        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4462        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4463    
4464        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4465    
4466        @param arg0: first argument of rank 2 or 4
4467        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4468        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4469        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4470        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4471        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4472        """
4473        sh0=getShape(arg0)
4474        sh1=getShape(arg1)
4475        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4476           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4477        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4478           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4479        else:
4480            raise ValueError,"first argument must have rank 2 or 4"
4481    
4482    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4483        """
4484        generalized tensor product of transposed of arg0 and arg1:
4485    
4486        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4487    
4488        where
4489    
4490            - s runs through arg0.Shape[axis_offset:]
4491            - r runs trough arg0.Shape[:axis_offset]
4492            - t runs through arg1.Shape[axis_offset:]
4493    
4494        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4495        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4496    
4497        @param arg0: first argument
4498        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4499        @param arg1: second argument
4500        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4501        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4502        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4503        """
4504        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4505        arg0,arg1=matchType(arg0,arg1)
4506        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4507        if isinstance(arg0,numarray.NumArray):
4508           if isinstance(arg1,Symbol):
4509               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4510           else:
4511               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4512                   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)
4513               arg0_c=arg0.copy()
4514               arg1_c=arg1.copy()
4515               sh0,sh1=arg0.shape,arg1.shape
4516               d0,d1,d01=1,1,1
4517               for i in sh0[axis_offset:]: d0*=i
4518               for i in sh1[axis_offset:]: d1*=i
4519               for i in sh0[:axis_offset]: d01*=i
4520               arg0_c.resize((d01,d0))
4521               arg1_c.resize((d01,d1))
4522               out=numarray.zeros((d0,d1),numarray.Float64)
4523               for i0 in range(d0):
4524                        for i1 in range(d1):
4525                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4526               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4527               return out
4528        elif isinstance(arg0,escript.Data):
4529           if isinstance(arg1,Symbol):
4530               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4531           else:
4532               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4533        else:      
4534           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4535                    
4536    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4537       """
4538       Symbol representing the general tensor product of the transposed of arg0 and arg1
4539       """
4540       def __init__(self,arg0,arg1,axis_offset=0):
4541           """
4542           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4543    
4544           @param arg0: first argument
4545           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4546           @param arg1: second argument
4547           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4548           @raise ValueError: inconsistent dimensions of arguments.
4549           @note: if both arguments have a spatial dimension, they must equal.
4550           """
4551           sh_arg0=getShape(arg0)
4552           sh_arg1=getShape(arg1)
4553           sh01=sh_arg0[:axis_offset]
4554           sh10=sh_arg1[:axis_offset]
4555           sh0=sh_arg0[axis_offset:]
4556           sh1=sh_arg1[axis_offset:]
4557           if not sh01==sh10:
4558               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)
4559           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4560    
4561       def getMyCode(self,argstrs,format="escript"):
4562          """
4563          returns a program code that can be used to evaluate the symbol.
4564    
4565          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4566          @type argstrs: C{list} of length 2 of C{str}.
4567          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4568          @type format: C{str}
4569          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4570          @rtype: C{str}
4571          @raise NotImplementedError: if the requested format is not available
4572          """
4573          if format=="escript" or format=="str" or format=="text":
4574             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4575          else:
4576             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4577    
4578       def substitute(self,argvals):
4579          """
4580          assigns new values to symbols in the definition of the symbol.
4581          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4582    
4583          @param argvals: new values assigned to symbols
4584          @type argvals: C{dict} with keywords of type L{Symbol}.
4585          @return: result of the substitution process. Operations are executed as much as possible.
4586          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4587          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4588          """
4589          if argvals.has_key(self):
4590             arg=argvals[self]
4591             if self.isAppropriateValue(arg):
4592                return arg
4593             else:
4594                raise TypeError,"%s: new value is not appropriate."%str(self)
4595          else:
4596             args=self.getSubstitutedArguments(argvals)
4597             return generalTransposedTensorProduct(args[0],args[1],args[2])
4598    
4599    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4600        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4601        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4602    
4603    def matrix_transposed_mult(arg0,arg1):
4604        """
4605        matrix-transposed(matrix) product of the two argument:
4606    
4607        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4608    
4609        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4610    
4611        The last dimensions of arg0 and arg1 must match.
4612    
4613        @param arg0: first argument of rank 2
4614        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4615        @param arg1: second argument of rank 2
4616        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4617        @return: the product of arg0 and the transposed of arg1 at each data point
4618        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4619        @raise ValueError: if the shapes of the arguments are not appropriate
4620        """
4621        sh0=getShape(arg0)
4622        sh1=getShape(arg1)
4623        if not len(sh0)==2 :
4624            raise ValueError,"first argument must have rank 2"
4625        if not len(sh1)==2 and not len(sh1)==1:
4626            raise ValueError,"second argument must have rank 1 or 2"
4627        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4628    
4629    def tensor_transposed_mult(arg0,arg1):
4630        """
4631        the tensor product of the first and the transpose of the second argument
4632        
4633        for arg0 of rank 2 this is
4634    
4635        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4636    
4637        and for arg0 of rank 4 this is
4638    
4639        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4640                  
4641        or
4642    
4643        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4644    
4645        In the first case the the second dimension of arg0 and arg1 must match and  
4646        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4647    
4648        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4649    
4650        @param arg0: first argument of rank 2 or 4
4651        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4652        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4653        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4654        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4655        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4656        """
4657        sh0=getShape(arg0)
4658        sh1=getShape(arg1)
4659        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4660           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4661        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4662           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4663        else:
4664            raise ValueError,"first argument must have rank 2 or 4"
4665    
4666    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4667        """
4668        generalized tensor product of transposed of arg0 and arg1:
4669    
4670        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4671    
4672        where
4673    
4674            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4675            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4676            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4677    
4678        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4679        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4680    
4681        @param arg0: first argument
4682        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4683        @param arg1: second argument
4684        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4685        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4686        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4687        """
4688        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4689        arg0,arg1=matchType(arg0,arg1)
4690        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4691        if isinstance(arg0,numarray.NumArray):
4692           if isinstance(arg1,Symbol):
4693               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4694           else:
4695               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4696                   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)
4697               arg0_c=arg0.copy()
4698               arg1_c=arg1.copy()
4699               sh0,sh1=arg0.shape,arg1.shape
4700               d0,d1,d01=1,1,1
4701               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4702               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4703               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4704               arg0_c.resize((d0,d01))
4705               arg1_c.resize((d1,d01))
4706               out=numarray.zeros((d0,d1),numarray.Float64)
4707               for i0 in range(d0):
4708                        for i1 in range(d1):
4709                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4710               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4711               return out
4712        elif isinstance(arg0,escript.Data):
4713           if isinstance(arg1,Symbol):
4714               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4715           else:
4716               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4717        else:      
4718           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4719                    
4720    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4721       """
4722       Symbol representing the general tensor product of arg0 and the transpose of arg1
4723       """
4724       def __init__(self,arg0,arg1,axis_offset=0):
4725           """
4726           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4727    
4728           @param arg0: first argument
4729           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4730           @param arg1: second argument
4731           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4732           @raise ValueError: inconsistent dimensions of arguments.
4733           @note: if both arguments have a spatial dimension, they must equal.
4734           """
4735           sh_arg0=getShape(arg0)
4736           sh_arg1=getShape(arg1)
4737           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4738           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4739           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4740           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4741           if not sh01==sh10:
4742               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)
4743           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4744    
4745       def getMyCode(self,argstrs,format="escript"):
4746          """
4747          returns a program code that can be used to evaluate the symbol.
4748    
4749          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4750          @type argstrs: C{list} of length 2 of C{str}.
4751          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4752          @type format: C{str}
4753          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4754          @rtype: C{str}
4755          @raise NotImplementedError: if the requested format is not available
4756          """
4757          if format=="escript" or format=="str" or format=="text":
4758             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4759          else:
4760             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4761    
4762       def substitute(self,argvals):
4763          """
4764          assigns new values to symbols in the definition of the symbol.
4765          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4766    
4767          @param argvals: new values assigned to symbols
4768          @type argvals: C{dict} with keywords of type L{Symbol}.
4769          @return: result of the substitution process. Operations are executed as much as possible.
4770          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4771          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4772          """
4773          if argvals.has_key(self):
4774             arg=argvals[self]
4775             if self.isAppropriateValue(arg):
4776                return arg
4777             else:
4778                raise TypeError,"%s: new value is not appropriate."%str(self)
4779          else:
4780             args=self.getSubstitutedArguments(argvals)
4781             return generalTensorTransposedProduct(args[0],args[1],args[2])
4782    
4783    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4784        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4785        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4786    
4787  #=========================================================  #=========================================================
4788  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency

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

  ViewVC Help
Powered by ViewVC 1.1.26