/[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 775 by ksteube, Mon Jul 10 04:00:08 2006 UTC
# Line 2966  def trace(arg,axis_offset=0): Line 2966  def trace(arg,axis_offset=0):
2966     else:     else:
2967        raise TypeError,"trace: Unknown argument type."        raise TypeError,"trace: Unknown argument type."
2968    
2969  def escript_trace(arg,axis_offset): # this should be escript._trace  def escript_trace(arg,axis_offset):
2970        "arg si a Data objects!!!"        "arg is a Data object"
2971        if arg.getRank()<2:        if arg.getRank()<2:
2972          raise ValueError,"escript_trace: rank of argument must be greater than 1"          raise ValueError,"escript_trace: rank of argument must be greater than 1"
2973        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
# Line 2975  def escript_trace(arg,axis_offset): # th Line 2975  def escript_trace(arg,axis_offset): # th
2975        s=list(arg.getShape())                s=list(arg.getShape())        
2976        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
2977          raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2978        out=escript.Data(0.,tuple(s[0:axis_offset]+s[axis_offset+2:]),arg.getFunctionSpace())        out=arg._matrixtrace(axis_offset)
       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]  
2979        return out        return out
2980  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
2981     """     """
# Line 3111  def transpose(arg,axis_offset=None): Line 3083  def transpose(arg,axis_offset=None):
3083     else:     else:
3084        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"transpose: Unknown argument type."
3085    
3086  def escript_transpose(arg,axis_offset): # this should be escript._transpose  def escript_transpose(arg,axis_offset):
3087        "arg si a Data objects!!!"        "arg is a Data object"
3088        r=arg.getRank()        r=arg.getRank()
3089        if axis_offset<0 or axis_offset>r:        if axis_offset<0 or axis_offset>r:
3090          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3091        s=arg.getShape()        out=arg._transpose(axis_offset)
       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.  
3092        return out        return out
3093  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3094     """     """
# Line 3284  def symmetric(arg): Line 3199  def symmetric(arg):
3199      else:      else:
3200        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3201    
3202  def escript_symmetric(arg): # this should be implemented in c++  def escript_symmetric(arg):
3203        if arg.getRank()==2:        if arg.getRank()==2:
3204          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3205             raise ValueError,"escript_symmetric: argument must be square."             raise ValueError,"escript_symmetric: argument must be square."
3206          out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())          out=arg._symmetric()
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               out[i0,i1]=(arg[i0,i1]+arg[i1,i0])/2.  
3207        elif arg.getRank()==4:        elif arg.getRank()==4:
3208          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]):
3209             raise ValueError,"escript_symmetric: argument must be square."             raise ValueError,"escript_symmetric: argument must be square."
3210          out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())          out=arg._symmetric()
         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.  
3211        else:        else:
3212          raise ValueError,"escript_symmetric: rank 2 or 4 is required."          raise ValueError,"escript_symmetric: rank 2 or 4 is required."
3213        return out        return out
# Line 3347  def escript_nonsymmetric(arg): # this sh Line 3254  def escript_nonsymmetric(arg): # this sh
3254        if arg.getRank()==2:        if arg.getRank()==2:
3255          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3256             raise ValueError,"escript_nonsymmetric: argument must be square."             raise ValueError,"escript_nonsymmetric: argument must be square."
3257          out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())          out=arg._nonsymmetric()
         for i0 in range(arg.getShape()[0]):  
            for i1 in range(arg.getShape()[1]):  
               out[i0,i1]=(arg[i0,i1]-arg[i1,i0])/2.  
3258        elif arg.getRank()==4:        elif arg.getRank()==4:
3259          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]):
3260             raise ValueError,"escript_nonsymmetric: argument must be square."             raise ValueError,"escript_nonsymmetric: argument must be square."
3261          out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())          out=arg._nonsymmetric()
         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.  
3262        else:        else:
3263          raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."          raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."
3264        return out        return out

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

  ViewVC Help
Powered by ViewVC 1.1.26