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

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