/[escript]/trunk-mpi-branch/escript/py_src/util.py
ViewVC logotype

Diff of /trunk-mpi-branch/escript/py_src/util.py

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

revision 574 by gross, Thu Mar 2 06:31:24 2006 UTC revision 580 by gross, Wed Mar 8 05:45:51 2006 UTC
# Line 3552  def eigenvalues(arg): Line 3552  def eigenvalues(arg):
3552      @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.      @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3553      """      """
3554      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3555        if not arg.rank==2:        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3556          raise ValueError,"eigenvalues: argument must have rank 2"        out.sort()
       s=arg.shape        
       if not s[0] == s[1]:  
         raise ValueError,"eigenvalues: argument must be a square matrix."  
       if s[0]==1:  
           out=arg[0]  
       elif s[0]==2:  
           A11=arg[0,0]  
           A12=arg[0,1]  
           A22=arg[1,1]  
           trA=(A11+A22)/2.  
           A11-=trA  
           A22-=trA  
           s=sqrt(A12**2-A11*A22)  
           out=trA+numarray.array([-s,s],type=numarray.Float64)  
       elif s[0]==3:  
           A11=arg[0,0]  
           A12=arg[0,1]  
           A22=arg[1,1]  
           A13=arg[0,2]  
           A23=arg[1,2]  
           A33=arg[2,2]  
           trA=(A11+A22+A33)/3.  
           A11-=trA  
           A22-=trA  
           A33-=trA  
           A13_2=A13**2  
           A23_2=A23**2  
           A12_2=A12**2  
           p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.  
           q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13  
           sq_p=sqrt(p/3.)  
           alpha_3=acos(clip(-q*sq_p**(-3.)/2.,-1.,1.))/3.  
           sq_p*=2.  
           out=trA+sq_p*numarray.array([-cos(alpha_3-numarray.pi/3.),-cos(alpha_3+numarray.pi/3.),cos(alpha_3)],type=numarray.Float64)  
       else:  
          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)  
          out.sort()  
3557        return out        return out
3558      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3559        return escript_eigenvalues(arg)        return arg._eigenvalues()
3560      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3561        if not arg.getRank()==2:        if not arg.getRank()==2:
3562          raise ValueError,"eigenvalues: argument must have rank 2"          raise ValueError,"eigenvalues: argument must have rank 2"
# Line 3642  def eigenvalues(arg): Line 3605  def eigenvalues(arg):
3605        return float(arg)        return float(arg)
3606      else:      else:
3607        raise TypeError,"eigenvalues: Unknown argument type."        raise TypeError,"eigenvalues: Unknown argument type."
   
 def escript_eigenvalues(arg): # this should be implemented in C++ arg and LAPACK is data object  
       if not arg.getRank()==2:  
         raise ValueError,"eigenvalues: argument must have rank 2"  
       s=arg.getShape()        
       if not s[0] == s[1]:  
         raise ValueError,"eigenvalues: argument must be a square matrix."  
       if s[0]==1:  
           return arg[0]  
       elif s[0]==2:  
           A11=arg[0,0]  
           A12=arg[0,1]  
           A22=arg[1,1]  
           trA=(A11+A22)/2.  
           A11-=trA  
           A22-=trA  
           s=sqrt(A12**2-A11*A22)  
           return trA+s*numarray.array([-1.,1.],type=numarray.Float64)  
       elif s[0]==3:  
           A11=arg[0,0]  
           A12=arg[0,1]  
           A22=arg[1,1]  
           A13=arg[0,2]  
           A23=arg[1,2]  
           A33=arg[2,2]  
           trA=(A11+A22+A33)/3.  
           A11-=trA  
           A22-=trA  
           A33-=trA  
           A13_2=A13**2  
           A23_2=A23**2  
           A12_2=A12**2  
           p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.  
           q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13  
           sq_p=sqrt(p/3.)  
           alpha_3=acos(clip(-q*sq_p**(-3.)/2.,-1.,1.))/3.  
           sq_p*=2.  
           f=escript.Data(0.,(3,),arg.getFunctionSpace())  
           f[0]=-cos(alpha_3-numarray.pi/3.)  
           f[1]=-cos(alpha_3+numarray.pi/3.)  
           f[2]=cos(alpha_3)  
           return trA+sq_p*f  
       else:  
          raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."  
3608  #=======================================================  #=======================================================
3609  #  Binary operations:  #  Binary operations:
3610  #=======================================================  #=======================================================

Legend:
Removed from v.574  
changed lines
  Added in v.580

  ViewVC Help
Powered by ViewVC 1.1.26