/[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 527 by gross, Tue Feb 14 02:25:02 2006 UTC revision 528 by gross, Wed Feb 15 02:20:33 2006 UTC
# Line 3282  def inverse(arg): Line 3282  def inverse(arg):
3282      """      """
3283      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3284    
3285      @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.
3286      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3287      @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])      @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3288      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3289        @remark: for L{escript.Data} objects the dimension is restricted to 3.
3290      """      """
3291      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3292        return numarray.linear_algebra.inverse(arg)        return numarray.linear_algebra.inverse(arg)
# Line 3420  class Inverse_Symbol(DependendSymbol): Line 3421  class Inverse_Symbol(DependendSymbol):
3421           return identity(self.getShape())           return identity(self.getShape())
3422        else:        else:
3423           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
3424    
3425    def eigenvalues(arg):
3426        """
3427        returns the eigenvalues of the square matrix arg.
3428    
3429        @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).
3431        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol},L {float} or L{int}
3432        @return: the list of the eigenvalues in increasing order.
3433        @rtype: C{list} of L{float}, 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.
3435        """
3436        if isinstance(arg,numarray.NumArray):
3437          if not arg.rank==2:
3438            raise ValueError,"eigenvalues: argument must have rank 2"
3439          s=arg.shape      
3440          if not s[0] == s[1]:
3441            raise ValueError,"eigenvalues: argument must be a square matrix."
3442          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3443          out.sort()
3444          return numarray.array2list(out)
3445        elif isinstance(arg,escript.Data) or isinstance(arg,Symbol):
3446          if not arg.getRank()==2:
3447            raise ValueError,"eigenvalues: argument must have rank 2"
3448          s=arg.getShape()      
3449          if not s[0] == s[1]:
3450            raise ValueError,"eigenvalues: argument must be a square matrix."
3451          if s[0]==1:
3452              return [arg[0,0]]
3453          elif s[0]==2:
3454              A11=arg[0,0]
3455              A12=arg[0,1]
3456              A22=arg[1,1]
3457              trA=(A11+A22)/2.
3458              A11-=trA
3459              A22-=trA
3460              s=sqrt(A12**2-A11*A22)
3461              return [trA-s,trA+s]
3462          elif s[0]==3:
3463              A11=arg[0,0]
3464              A12=arg[0,1]
3465              A22=arg[1,1]
3466              A13=arg[0,2]
3467              A23=arg[1,2]
3468              A33=arg[2,2]
3469              trA=(A11+A22+A22)/3.
3470              A11-=trA
3471              A22-=trA
3472              A33-=trA
3473              A13_2=A13**2
3474              A23_2=A23**2
3475              A12_2=A12**2
3476              p=A13_2+A_23_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
3478              sq_p=sqrt(p/3.)
3479              alpha_3=acos(-q/sq_p**(1./3.)/2.)/3.
3480              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.)]
3482          else:
3483             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3484        elif isinstance(arg,float):
3485          return [arg]
3486        elif isinstance(arg,int):
3487          return [float(arg)]
3488        else:
3489          raise TypeError,"eigenvalues: Unknown argument type."
3490  #=======================================================  #=======================================================
3491  #  Binary operations:  #  Binary operations:
3492  #=======================================================  #=======================================================
Legend:
Removed from v.527  
changed lines
  Added in v.528
  ViewVC Help
Powered by ViewVC 1.1.26