/[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 720 by gross, Thu Apr 27 10:16:05 2006 UTC revision 795 by ksteube, Mon Jul 31 01:23:58 2006 UTC
# Line 2939  def trace(arg,axis_offset=0): Line 2939  def trace(arg,axis_offset=0):
2939     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
2940        sh=arg.shape        sh=arg.shape
2941        if len(sh)<2:        if len(sh)<2:
2942          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2943        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
2944          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
2945        s1=1        s1=1
2946        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
2947        s2=1        s2=1
2948        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2949        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
2950          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)
2951        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))
2952        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
2953        for i1 in range(s1):        for i1 in range(s1):
# Line 2956  def trace(arg,axis_offset=0): Line 2956  def trace(arg,axis_offset=0):
2956        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2957        return out        return out
2958     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
2959        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
2960            raise ValueError,"rank of argument must be greater than 1"
2961          if axis_offset<0 or axis_offset>arg.getRank()-2:
2962            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2963          s=list(arg.getShape())        
2964          if not s[axis_offset] == s[axis_offset+1]:
2965            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2966          return arg._matrixtrace(axis_offset)
2967     elif isinstance(arg,float):     elif isinstance(arg,float):
2968        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
2969     elif isinstance(arg,int):     elif isinstance(arg,int):
2970        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
2971     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
2972        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
2973     else:     else:
2974        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
2975    
 def escript_trace(arg,axis_offset): # this should be escript._trace  
       "arg si a Data objects!!!"  
       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=escript.Data(0.,tuple(s[0:axis_offset]+s[axis_offset+2:]),arg.getFunctionSpace())  
       if arg.getRank()==2:  
          for i0 in range(s[0]):  
             out+=arg[i0,i0]  
       elif arg.getRank()==3:  
          if axis_offset==0:  
             for i0 in range(s[0]):  
                   for i2 in range(s[2]):  
                          out[i2]+=arg[i0,i0,i2]  
          elif axis_offset==1:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                          out[i0]+=arg[i0,i1,i1]  
       elif arg.getRank()==4:  
          if axis_offset==0:  
             for i0 in range(s[0]):  
                   for i2 in range(s[2]):  
                      for i3 in range(s[3]):  
                          out[i2,i3]+=arg[i0,i0,i2,i3]  
          elif axis_offset==1:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                      for i3 in range(s[3]):  
                          out[i0,i3]+=arg[i0,i1,i1,i3]  
          elif axis_offset==2:  
             for i0 in range(s[0]):  
                for i1 in range(s[1]):  
                   for i2 in range(s[2]):  
                          out[i0,i1]+=arg[i0,i1,i2,i2]  
       return out  
2976  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
2977     """     """
2978     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 3019  class Trace_Symbol(DependendSymbol): Line 2987  class Trace_Symbol(DependendSymbol):
2987        @type axis_offset: C{int}        @type axis_offset: C{int}
2988        """        """
2989        if arg.getRank()<2:        if arg.getRank()<2:
2990          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2991        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
2992          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
2993        s=list(arg.getShape())                s=list(arg.getShape())        
2994        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
2995          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)
2996        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())
2997    
2998     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3095  def transpose(arg,axis_offset=None): Line 3063  def transpose(arg,axis_offset=None):
3063        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3064        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))
3065     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3066        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3067        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3068          if axis_offset<0 or axis_offset>r:
3069            raise ValueError,"axis_offset must be between 0 and %s"%r
3070          return arg._transpose(axis_offset)
3071     elif isinstance(arg,float):     elif isinstance(arg,float):
3072        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3073          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3074        return arg        return arg
3075     elif isinstance(arg,int):     elif isinstance(arg,int):
3076        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3077          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3078        return float(arg)        return float(arg)
3079     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3080        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3081        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3082     else:     else:
3083        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3084    
 def escript_transpose(arg,axis_offset): # this should be escript._transpose  
       "arg si a Data objects!!!"  
       r=arg.getRank()  
       if axis_offset<0 or axis_offset>r:  
         raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r  
       s=arg.getShape()  
       s_out=s[axis_offset:]+s[:axis_offset]  
       out=escript.Data(0.,s_out,arg.getFunctionSpace())  
       if r==4:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i3,i0,i1,i2]  
          elif axis_offset==2:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i2,i3,i0,i1]  
          elif axis_offset==3:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i1,i2,i3,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                      for i3 in range(s_out[3]):  
                          out[i0,i1,i2,i3]=arg[i0,i1,i2,i3]  
       elif r==3:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i2,i0,i1]  
          elif axis_offset==2:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i1,i2,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                   for i2 in range(s_out[2]):  
                          out[i0,i1,i2]=arg[i0,i1,i2]  
       elif r==2:  
          if axis_offset==1:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                          out[i0,i1]=arg[i1,i0]  
          else:  
             for i0 in range(s_out[0]):  
                for i1 in range(s_out[1]):  
                          out[i0,i1]=arg[i0,i1]  
       elif r==1:  
           for i0 in range(s_out[0]):  
                out[i0]=arg[i0]  
       elif r==0:  
              out=arg+0.  
       return out  
3085  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3086     """     """
3087     L{Symbol} representing the result of the transpose function     L{Symbol} representing the result of the transpose function
# Line 3191  class Transpose_Symbol(DependendSymbol): Line 3098  class Transpose_Symbol(DependendSymbol):
3098        """        """
3099        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3100        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3101          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%r
3102        s=arg.getShape()        s=arg.getShape()
3103        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())
3104    
# Line 3258  def symmetric(arg): Line 3165  def symmetric(arg):
3165      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3166        if arg.rank==2:        if arg.rank==2:
3167          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3168             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3169        elif arg.rank==4:        elif arg.rank==4:
3170          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]):
3171             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3172        else:        else:
3173          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3174        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3175      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3176        return escript_symmetric(arg)        if arg.getRank()==2:
3177            if not (arg.getShape()[0]==arg.getShape()[1]):
3178               raise ValueError,"argument must be square."
3179            return arg._symmetric()
3180          elif arg.getRank()==4:
3181            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3182               raise ValueError,"argument must be square."
3183            return arg._symmetric()
3184          else:
3185            raise ValueError,"rank 2 or 4 is required."
3186      elif isinstance(arg,float):      elif isinstance(arg,float):
3187        return arg        return arg
3188      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3274  def symmetric(arg): Line 3190  def symmetric(arg):
3190      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3191        if arg.getRank()==2:        if arg.getRank()==2:
3192          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3193             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3194        elif arg.getRank()==4:        elif arg.getRank()==4:
3195          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]):
3196             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3197        else:        else:
3198          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3199        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3200      else:      else:
3201        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3202    
 def escript_symmetric(arg): # this should be implemented in c++  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               out[i0,i1]=(arg[i0,i1]+arg[i1,i0])/2.  
       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=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               for i2 in range(arg.getShape()[2]):  
                  for i3 in range(arg.getShape()[3]):  
                      out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]+arg[i2,i3,i0,i1])/2.  
       else:  
         raise ValueError,"escript_symmetric: rank 2 or 4 is required."  
       return out  
   
3203  def nonsymmetric(arg):  def nonsymmetric(arg):
3204      """      """
3205      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 3325  def nonsymmetric(arg): Line 3220  def nonsymmetric(arg):
3220          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3221        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3222      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3223        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3224            if not (arg.getShape()[0]==arg.getShape()[1]):
3225               raise ValueError,"argument must be square."
3226            return arg._nonsymmetric()
3227          elif arg.getRank()==4:
3228            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3229               raise ValueError,"argument must be square."
3230            return arg._nonsymmetric()
3231          else:
3232            raise ValueError,"rank 2 or 4 is required."
3233      elif isinstance(arg,float):      elif isinstance(arg,float):
3234        return arg        return arg
3235      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3343  def nonsymmetric(arg): Line 3247  def nonsymmetric(arg):
3247      else:      else:
3248        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3249    
 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=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               out[i0,i1]=(arg[i0,i1]-arg[i1,i0])/2.  
       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=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())  
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               for i2 in range(arg.getShape()[2]):  
                  for i3 in range(arg.getShape()[3]):  
                      out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]-arg[i2,i3,i0,i1])/2.  
       else:  
         raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."  
       return out  
   
   
3250  def inverse(arg):  def inverse(arg):
3251      """      """
3252      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3253    
3254      @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.
3255      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3256      @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])
3257      @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
3258      @note: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3259      """      """
# Line 3508  class Inverse_Symbol(DependendSymbol): Line 3390  class Inverse_Symbol(DependendSymbol):
3390        if arg==self:        if arg==self:
3391           return identity(self.getShape())           return identity(self.getShape())
3392        else:        else:
3393           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3394    
3395  def eigenvalues(arg):  def eigenvalues(arg):
3396      """      """
# Line 4096  def inner(arg0,arg1): Line 3978  def inner(arg0,arg1):
3978          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
3979      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
3980    
3981    def outer(arg0,arg1):
3982        """
3983        the outer product of the two argument:
3984    
3985        out[t,s]=arg0[t]*arg1[s]
3986    
3987        where
3988    
3989            - s runs through arg0.Shape
3990            - t runs through arg1.Shape
3991    
3992        @param arg0: first argument
3993        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
3994        @param arg1: second argument
3995        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
3996        @return: the outer product of arg0 and arg1 at each data point
3997        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3998        """
3999        return generalTensorProduct(arg0,arg1,axis_offset=0)
4000    
4001  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4002      """      """
4003        see L{matrix_mult}
4004        """
4005        return matrix_mult(arg0,arg1)
4006    
4007    def matrix_mult(arg0,arg1):
4008        """
4009      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4010    
4011      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
# Line 4106  def matrixmult(arg0,arg1): Line 4014  def matrixmult(arg0,arg1):
4014    
4015      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4016    
4017      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.
4018    
4019      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4020      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4124  def matrixmult(arg0,arg1): Line 4032  def matrixmult(arg0,arg1):
4032          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4033      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4034    
4035  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4036      """      """
4037      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  
4038      """      """
4039      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
   
4040    
4041  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4042      """      """
4043      the tensor product of the two argument:      the tensor product of the two argument:
4044            
# Line 4170  def tensormult(arg0,arg1): Line 4063  def tensormult(arg0,arg1):
4063    
4064      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]
4065    
4066      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  
4067      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.
4068    
4069      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4070      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4187  def tensormult(arg0,arg1): Line 4080  def tensormult(arg0,arg1):
4080      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):
4081         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4082      else:      else:
4083          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4084    
4085  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4086      """      """
# Line 4201  def generalTensorProduct(arg0,arg1,axis_ Line 4094  def generalTensorProduct(arg0,arg1,axis_
4094          - r runs trough arg0.Shape[:axis_offset]          - r runs trough arg0.Shape[:axis_offset]
4095          - t runs through arg1.Shape[axis_offset:]          - t runs through arg1.Shape[axis_offset:]
4096    
     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.  
   
4097      @param arg0: first argument      @param arg0: first argument
4098      @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}
4099      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4100      @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}
4101      @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.
4102      @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 4219  def generalTensorProduct(arg0,arg1,axis_ Line 4109  def generalTensorProduct(arg0,arg1,axis_
4109             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4110         else:         else:
4111             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4112                 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)
4113             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4114             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4115             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4245  def generalTensorProduct(arg0,arg1,axis_ Line 4135  def generalTensorProduct(arg0,arg1,axis_
4135                                    
4136  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4137     """     """
4138     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4139     """     """
4140     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4141         """         """
4142         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4143    
4144         @param arg0: numerator         @param arg0: first argument
4145         @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}.
4146         @param arg1: denominator         @param arg1: second argument
4147         @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}.
4148         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4149         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4150         """         """
4151         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
# Line 4306  class GeneralTensorProduct_Symbol(Depend Line 4196  class GeneralTensorProduct_Symbol(Depend
4196           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4197           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4198    
4199    def escript_generalTensorProductNew(arg0,arg1,axis_offset):
4200        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4201        # calculate the return shape:
4202        shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4203        shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4204        shape10=arg1.getShape()[:axis_offset]
4205        shape1=arg1.getShape()[axis_offset:]
4206        if not shape01==shape10:
4207            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)
4208        # Figure out which functionspace to use (look at where operator+ is defined maybe in BinaryOp.h to get the logic for this)
4209        # fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4210        out=GeneralTensorProduct(arg0, arg1, axis_offset)
4211        return out
4212    
4213  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct
4214      "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!!!"
4215      # calculate the return shape:      # calculate the return shape:
# Line 4348  def escript_generalTensorProduct(arg0,ar Line 4252  def escript_generalTensorProduct(arg0,ar
4252           out.__setitem__(tuple(i0+i1),s)           out.__setitem__(tuple(i0+i1),s)
4253      return out      return out
4254    
4255    
4256    def transposed_matrix_mult(arg0,arg1):
4257        """
4258        transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4259    
4260        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4261    
4262        or
4263    
4264        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4265    
4266        The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4267    
4268        The first dimension of arg0 and arg1 must match.
4269    
4270        @param arg0: first argument of rank 2
4271        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4272        @param arg1: second argument of at least rank 1
4273        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4274        @return: the product of the transposed of arg0 and arg1 at each data point
4275        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4276        @raise ValueError: if the shapes of the arguments are not appropriate
4277        """
4278        sh0=pokeShape(arg0)
4279        sh1=pokeShape(arg1)
4280        if not len(sh0)==2 :
4281            raise ValueError,"first argument must have rank 2"
4282        if not len(sh1)==2 and not len(sh1)==1:
4283            raise ValueError,"second argument must have rank 1 or 2"
4284        return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4285    
4286    def transposed_tensor_mult(arg0,arg1):
4287        """
4288        the tensor product of the transposed of the first and the second argument
4289        
4290        for arg0 of rank 2 this is
4291    
4292        out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4293    
4294        or
4295    
4296        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4297    
4298      
4299        and for arg0 of rank 4 this is
4300    
4301        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4302                  
4303        or
4304    
4305        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4306    
4307        or
4308    
4309        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4310    
4311        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4312        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4313    
4314        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4315    
4316        @param arg0: first argument of rank 2 or 4
4317        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4318        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4319        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4320        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4321        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4322        """
4323        sh0=pokeShape(arg0)
4324        sh1=pokeShape(arg1)
4325        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4326           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4327        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4328           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4329        else:
4330            raise ValueError,"first argument must have rank 2 or 4"
4331    
4332    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4333        """
4334        generalized tensor product of transposed of arg0 and arg1:
4335    
4336        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4337    
4338        where
4339    
4340            - s runs through arg0.Shape[axis_offset:]
4341            - r runs trough arg0.Shape[:axis_offset]
4342            - t runs through arg1.Shape[axis_offset:]
4343    
4344        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4345        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4346    
4347        @param arg0: first argument
4348        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4349        @param arg1: second argument
4350        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4351        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4352        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4353        """
4354        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4355        arg0,arg1=matchType(arg0,arg1)
4356        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4357        if isinstance(arg0,numarray.NumArray):
4358           if isinstance(arg1,Symbol):
4359               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4360           else:
4361               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4362                   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)
4363               arg0_c=arg0.copy()
4364               arg1_c=arg1.copy()
4365               sh0,sh1=arg0.shape,arg1.shape
4366               d0,d1,d01=1,1,1
4367               for i in sh0[axis_offset:]: d0*=i
4368               for i in sh1[axis_offset:]: d1*=i
4369               for i in sh0[:axis_offset]: d01*=i
4370               arg0_c.resize((d01,d0))
4371               arg1_c.resize((d01,d1))
4372               out=numarray.zeros((d0,d1),numarray.Float64)
4373               for i0 in range(d0):
4374                        for i1 in range(d1):
4375                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4376               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4377               return out
4378        elif isinstance(arg0,escript.Data):
4379           if isinstance(arg1,Symbol):
4380               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4381           else:
4382               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4383        else:      
4384           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4385                    
4386    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4387       """
4388       Symbol representing the general tensor product of the transposed of arg0 and arg1
4389       """
4390       def __init__(self,arg0,arg1,axis_offset=0):
4391           """
4392           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4393    
4394           @param arg0: first argument
4395           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4396           @param arg1: second argument
4397           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4398           @raise ValueError: inconsistent dimensions of arguments.
4399           @note: if both arguments have a spatial dimension, they must equal.
4400           """
4401           sh_arg0=pokeShape(arg0)
4402           sh_arg1=pokeShape(arg1)
4403           sh01=sh_arg0[:axis_offset]
4404           sh10=sh_arg1[:axis_offset]
4405           sh0=sh_arg0[axis_offset:]
4406           sh1=sh_arg1[axis_offset:]
4407           if not sh01==sh10:
4408               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)
4409           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4410    
4411       def getMyCode(self,argstrs,format="escript"):
4412          """
4413          returns a program code that can be used to evaluate the symbol.
4414    
4415          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4416          @type argstrs: C{list} of length 2 of C{str}.
4417          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4418          @type format: C{str}
4419          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4420          @rtype: C{str}
4421          @raise NotImplementedError: if the requested format is not available
4422          """
4423          if format=="escript" or format=="str" or format=="text":
4424             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4425          else:
4426             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4427    
4428       def substitute(self,argvals):
4429          """
4430          assigns new values to symbols in the definition of the symbol.
4431          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4432    
4433          @param argvals: new values assigned to symbols
4434          @type argvals: C{dict} with keywords of type L{Symbol}.
4435          @return: result of the substitution process. Operations are executed as much as possible.
4436          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4437          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4438          """
4439          if argvals.has_key(self):
4440             arg=argvals[self]
4441             if self.isAppropriateValue(arg):
4442                return arg
4443             else:
4444                raise TypeError,"%s: new value is not appropriate."%str(self)
4445          else:
4446             args=self.getSubstitutedArguments(argvals)
4447             return generalTransposedTensorProduct(args[0],args[1],args[2])
4448    
4449    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4450        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4451        # calculate the return shape:
4452        shape01=arg0.getShape()[:axis_offset]
4453        shape10=arg1.getShape()[:axis_offset]
4454        shape0=arg0.getShape()[axis_offset:]
4455        shape1=arg1.getShape()[axis_offset:]
4456        if not shape01==shape10:
4457            raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4458    
4459        # whatr function space should be used? (this here is not good!)
4460        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4461        # create return value:
4462        out=escript.Data(0.,tuple(shape0+shape1),fs)
4463        #
4464        s0=[[]]
4465        for k in shape0:
4466              s=[]
4467              for j in s0:
4468                    for i in range(k): s.append(j+[slice(i,i)])
4469              s0=s
4470        s1=[[]]
4471        for k in shape1:
4472              s=[]
4473              for j in s1:
4474                    for i in range(k): s.append(j+[slice(i,i)])
4475              s1=s
4476        s01=[[]]
4477        for k in shape01:
4478              s=[]
4479              for j in s01:
4480                    for i in range(k): s.append(j+[slice(i,i)])
4481              s01=s
4482    
4483        for i0 in s0:
4484           for i1 in s1:
4485             s=escript.Scalar(0.,fs)
4486             for i01 in s01:
4487                s+=arg0.__getitem__(tuple(i01+i0))*arg1.__getitem__(tuple(i01+i1))
4488             out.__setitem__(tuple(i0+i1),s)
4489        return out
4490    
4491    
4492    def matrix_transposed_mult(arg0,arg1):
4493        """
4494        matrix-transposed(matrix) product of the two argument:
4495    
4496        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4497    
4498        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4499    
4500        The last dimensions of arg0 and arg1 must match.
4501    
4502        @param arg0: first argument of rank 2
4503        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4504        @param arg1: second argument of rank 2
4505        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4506        @return: the product of arg0 and the transposed of arg1 at each data point
4507        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4508        @raise ValueError: if the shapes of the arguments are not appropriate
4509        """
4510        sh0=pokeShape(arg0)
4511        sh1=pokeShape(arg1)
4512        if not len(sh0)==2 :
4513            raise ValueError,"first argument must have rank 2"
4514        if not len(sh1)==2 and not len(sh1)==1:
4515            raise ValueError,"second argument must have rank 1 or 2"
4516        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4517    
4518    def tensor_transposed_mult(arg0,arg1):
4519        """
4520        the tensor product of the first and the transpose of the second argument
4521        
4522        for arg0 of rank 2 this is
4523    
4524        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4525    
4526        and for arg0 of rank 4 this is
4527    
4528        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4529                  
4530        or
4531    
4532        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4533    
4534        In the first case the the second dimension of arg0 and arg1 must match and  
4535        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4536    
4537        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4538    
4539        @param arg0: first argument of rank 2 or 4
4540        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4541        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4542        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4543        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4544        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4545        """
4546        sh0=pokeShape(arg0)
4547        sh1=pokeShape(arg1)
4548        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4549           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4550        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4551           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4552        else:
4553            raise ValueError,"first argument must have rank 2 or 4"
4554    
4555    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4556        """
4557        generalized tensor product of transposed of arg0 and arg1:
4558    
4559        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4560    
4561        where
4562    
4563            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4564            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4565            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4566    
4567        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4568        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4569    
4570        @param arg0: first argument
4571        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4572        @param arg1: second argument
4573        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4574        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4575        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4576        """
4577        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4578        arg0,arg1=matchType(arg0,arg1)
4579        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4580        if isinstance(arg0,numarray.NumArray):
4581           if isinstance(arg1,Symbol):
4582               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4583           else:
4584               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4585                   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)
4586               arg0_c=arg0.copy()
4587               arg1_c=arg1.copy()
4588               sh0,sh1=arg0.shape,arg1.shape
4589               d0,d1,d01=1,1,1
4590               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4591               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4592               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4593               arg0_c.resize((d0,d01))
4594               arg1_c.resize((d1,d01))
4595               out=numarray.zeros((d0,d1),numarray.Float64)
4596               for i0 in range(d0):
4597                        for i1 in range(d1):
4598                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4599               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4600               return out
4601        elif isinstance(arg0,escript.Data):
4602           if isinstance(arg1,Symbol):
4603               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4604           else:
4605               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4606        else:      
4607           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4608                    
4609    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4610       """
4611       Symbol representing the general tensor product of arg0 and the transpose of arg1
4612       """
4613       def __init__(self,arg0,arg1,axis_offset=0):
4614           """
4615           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4616    
4617           @param arg0: first argument
4618           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4619           @param arg1: second argument
4620           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4621           @raise ValueError: inconsistent dimensions of arguments.
4622           @note: if both arguments have a spatial dimension, they must equal.
4623           """
4624           sh_arg0=pokeShape(arg0)
4625           sh_arg1=pokeShape(arg1)
4626           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4627           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4628           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4629           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4630           if not sh01==sh10:
4631               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)
4632           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4633    
4634       def getMyCode(self,argstrs,format="escript"):
4635          """
4636          returns a program code that can be used to evaluate the symbol.
4637    
4638          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4639          @type argstrs: C{list} of length 2 of C{str}.
4640          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4641          @type format: C{str}
4642          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4643          @rtype: C{str}
4644          @raise NotImplementedError: if the requested format is not available
4645          """
4646          if format=="escript" or format=="str" or format=="text":
4647             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4648          else:
4649             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4650    
4651       def substitute(self,argvals):
4652          """
4653          assigns new values to symbols in the definition of the symbol.
4654          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4655    
4656          @param argvals: new values assigned to symbols
4657          @type argvals: C{dict} with keywords of type L{Symbol}.
4658          @return: result of the substitution process. Operations are executed as much as possible.
4659          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4660          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4661          """
4662          if argvals.has_key(self):
4663             arg=argvals[self]
4664             if self.isAppropriateValue(arg):
4665                return arg
4666             else:
4667                raise TypeError,"%s: new value is not appropriate."%str(self)
4668          else:
4669             args=self.getSubstitutedArguments(argvals)
4670             return generalTensorTransposedProduct(args[0],args[1],args[2])
4671    
4672    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4673        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4674        # calculate the return shape:
4675        shape01=arg0.getShape()[arg0.getRank()-axis_offset:]
4676        shape0=arg0.getShape()[:arg0.getRank()-axis_offset]
4677        shape10=arg1.getShape()[arg1.getRank()-axis_offset:]
4678        shape1=arg1.getShape()[:arg1.getRank()-axis_offset]
4679        if not shape01==shape10:
4680            raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
4681    
4682        # whatr function space should be used? (this here is not good!)
4683        fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()
4684        # create return value:
4685        out=escript.Data(0.,tuple(shape0+shape1),fs)
4686        #
4687        s0=[[]]
4688        for k in shape0:
4689              s=[]
4690              for j in s0:
4691                    for i in range(k): s.append(j+[slice(i,i)])
4692              s0=s
4693        s1=[[]]
4694        for k in shape1:
4695              s=[]
4696              for j in s1:
4697                    for i in range(k): s.append(j+[slice(i,i)])
4698              s1=s
4699        s01=[[]]
4700        for k in shape01:
4701              s=[]
4702              for j in s01:
4703                    for i in range(k): s.append(j+[slice(i,i)])
4704              s01=s
4705    
4706        for i0 in s0:
4707           for i1 in s1:
4708             s=escript.Scalar(0.,fs)
4709             for i01 in s01:
4710                s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i1+i01))
4711             out.__setitem__(tuple(i0+i1),s)
4712        return out
4713    
4714    
4715  #=========================================================  #=========================================================
4716  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency

Legend:
Removed from v.720  
changed lines
  Added in v.795

  ViewVC Help
Powered by ViewVC 1.1.26