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) |
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 |
#======================================================= |
#======================================================= |