/[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 528 by gross, Wed Feb 15 02:20:33 2006 UTC revision 529 by gross, Wed Feb 15 03:58:50 2006 UTC
# Line 3428  def eigenvalues(arg): Line 3428  def eigenvalues(arg):
3428    
3429      @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.
3430                  arg must be symmetric, ie. transpose(arg)==arg (this is not checked).                  arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3431      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol},L {float} or L{int}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3432      @return: the list of the eigenvalues in increasing order.      @return: the eigenvalues in increasing order.
3433      @rtype: C{list} of L{float}, L{escript.Data}, L{Symbol} depending on the input.      @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3434      @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.
3435      """      """
3436      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
# Line 3441  def eigenvalues(arg): Line 3441  def eigenvalues(arg):
3441          raise ValueError,"eigenvalues: argument must be a square matrix."          raise ValueError,"eigenvalues: argument must be a square matrix."
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 numarray.array2list(out)        return out
3445      elif isinstance(arg,escript.Data) or isinstance(arg,Symbol):      elif isinstance(arg,escript.Data) or isinstance(arg,Symbol):
3446        if not arg.getRank()==2:        if not arg.getRank()==2:
3447          raise ValueError,"eigenvalues: argument must have rank 2"          raise ValueError,"eigenvalues: argument must have rank 2"
# Line 3449  def eigenvalues(arg): Line 3449  def eigenvalues(arg):
3449        if not s[0] == s[1]:        if not s[0] == s[1]:
3450          raise ValueError,"eigenvalues: argument must be a square matrix."          raise ValueError,"eigenvalues: argument must be a square matrix."
3451        if s[0]==1:        if s[0]==1:
3452            return [arg[0,0]]            return arg[0]
3453        elif s[0]==2:        elif s[0]==2:
3454            A11=arg[0,0]            A11=arg[0,0]
3455            A12=arg[0,1]            A12=arg[0,1]
# Line 3458  def eigenvalues(arg): Line 3458  def eigenvalues(arg):
3458            A11-=trA            A11-=trA
3459            A22-=trA            A22-=trA
3460            s=sqrt(A12**2-A11*A22)            s=sqrt(A12**2-A11*A22)
3461            return [trA-s,trA+s]            return trA+s*numarray.array([-1.,1.])
3462        elif s[0]==3:        elif s[0]==3:
3463            A11=arg[0,0]            A11=arg[0,0]
3464            A12=arg[0,1]            A12=arg[0,1]
# Line 3473  def eigenvalues(arg): Line 3473  def eigenvalues(arg):
3473            A13_2=A13**2            A13_2=A13**2
3474            A23_2=A23**2            A23_2=A23**2
3475            A12_2=A12**2            A12_2=A12**2
3476            p=A13_2+A_23_2+A12_2+(A11**2+A22**2+A33**2)/2.            p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3477            q=A13_2*A22+A_23_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
3478            sq_p=sqrt(p/3.)            sq_p=sqrt(p/3.)
3479            alpha_3=acos(-q/sq_p**(1./3.)/2.)/3.            alpha_3=acos(-q/sq_p**(1./3.)/2.)/3.
3480            sq_p*=2.            sq_p*=2.
3481            return [trA+sq_p*cos(alpha_3),trA-sq_p*cos(alpha_3+numarray.pi/3.),trA-sq_p*cos(alpha_3-numarray.pi/3.)]            f=cos(alpha_3)               *numarray.array([1.,0.,0.]) \
3482               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.]) \
3483               -cos(alpha_3-numarray.pi/3.)*numarray.array([0.,0.,1.])
3484              return trA+sq_p*f
3485        else:        else:
3486           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."
3487      elif isinstance(arg,float):      elif isinstance(arg,float):
3488        return [arg]        return arg
3489      elif isinstance(arg,int):      elif isinstance(arg,int):
3490        return [float(arg)]        return float(arg)
3491      else:      else:
3492        raise TypeError,"eigenvalues: Unknown argument type."        raise TypeError,"eigenvalues: Unknown argument type."
3493  #=======================================================  #=======================================================
Legend:
Removed from v.528  
changed lines
  Added in v.529
  ViewVC Help
Powered by ViewVC 1.1.26