/[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 529 by gross, Wed Feb 15 03:58:50 2006 UTC revision 530 by gross, Wed Feb 15 07:11:00 2006 UTC
# Line 3442  def eigenvalues(arg): Line 3442  def eigenvalues(arg):
3442        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)        out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3443        out.sort()        out.sort()
3444        return out        return out
3445      elif isinstance(arg,escript.Data) or isinstance(arg,Symbol):      elif isinstance(arg,escript.Data):
3446          return escript_eigenvalues(arg)
3447        elif isinstance(arg,Symbol):
3448        if not arg.getRank()==2:        if not arg.getRank()==2:
3449          raise ValueError,"eigenvalues: argument must have rank 2"          raise ValueError,"eigenvalues: argument must have rank 2"
3450        s=arg.getShape()              s=arg.getShape()      
# Line 3466  def eigenvalues(arg): Line 3468  def eigenvalues(arg):
3468            A13=arg[0,2]            A13=arg[0,2]
3469            A23=arg[1,2]            A23=arg[1,2]
3470            A33=arg[2,2]            A33=arg[2,2]
3471            trA=(A11+A22+A22)/3.            trA=(A11+A22+A33)/3.
3472            A11-=trA            A11-=trA
3473            A22-=trA            A22-=trA
3474            A33-=trA            A33-=trA
# Line 3476  def eigenvalues(arg): Line 3478  def eigenvalues(arg):
3478            p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.            p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3479            q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13            q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3480            sq_p=sqrt(p/3.)            sq_p=sqrt(p/3.)
3481            alpha_3=acos(-q/sq_p**(1./3.)/2.)/3.            alpha_3=acos(-q*sq_p**(-3.)/2.)/3.
3482            sq_p*=2.            sq_p*=2.
3483            f=cos(alpha_3)               *numarray.array([1.,0.,0.]) \            f=cos(alpha_3)               *numarray.array([0.,0.,1.]) \
3484             -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.]) \             -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.]) \
3485             -cos(alpha_3-numarray.pi/3.)*numarray.array([0.,0.,1.])             -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.])
3486            return trA+sq_p*f            return trA+sq_p*f
3487        else:        else:
3488           raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."           raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
# Line 3490  def eigenvalues(arg): Line 3492  def eigenvalues(arg):
3492        return float(arg)        return float(arg)
3493      else:      else:
3494        raise TypeError,"eigenvalues: Unknown argument type."        raise TypeError,"eigenvalues: Unknown argument type."
3495    
3496    def escript_eigenvalues(arg): # this should be implemented in C++ arg and LAPACK is data object
3497          if not arg.getRank()==2:
3498            raise ValueError,"eigenvalues: argument must have rank 2"
3499          s=arg.getShape()      
3500          if not s[0] == s[1]:
3501            raise ValueError,"eigenvalues: argument must be a square matrix."
3502          if s[0]==1:
3503              return arg[0]
3504          elif s[0]==2:
3505              A11=arg[0,0]
3506              A12=arg[0,1]
3507              A22=arg[1,1]
3508              trA=(A11+A22)/2.
3509              A11-=trA
3510              A22-=trA
3511              s=sqrt(A12**2-A11*A22)
3512              return trA+s*numarray.array([-1.,1.])
3513          elif s[0]==3:
3514              A11=arg[0,0]
3515              A12=arg[0,1]
3516              A22=arg[1,1]
3517              A13=arg[0,2]
3518              A23=arg[1,2]
3519              A33=arg[2,2]
3520              trA=(A11+A22+A33)/3.
3521              A11-=trA
3522              A22-=trA
3523              A33-=trA
3524              A13_2=A13**2
3525              A23_2=A23**2
3526              A12_2=A12**2
3527              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3528              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3529              sq_p=sqrt(p/3.)
3530              alpha_3=acos(-q*sq_p**(-3.)/2.)/3.
3531              sq_p*=2.
3532              f=escript.Data(0.,(3,),arg.getFunctionSpace())
3533              f[0]=-cos(alpha_3-numarray.pi/3.)
3534              f[1]=-cos(alpha_3+numarray.pi/3.)
3535              f[2]=cos(alpha_3)
3536              return trA+sq_p*f
3537          else:
3538             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3539  #=======================================================  #=======================================================
3540  #  Binary operations:  #  Binary operations:
3541  #=======================================================  #=======================================================
Legend:
Removed from v.529  
changed lines
  Added in v.530
  ViewVC Help
Powered by ViewVC 1.1.26