comp.lang.idl-pvwave archive
Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
Home » Public Forums » archive » IDL inverse matrix problem
Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend  		Switch to threaded view of this topic Create a new topic Submit Reply
IDL inverse matrix problem [message #89975] Mon, 12 January 2015 09:56 Go to next message
amin farhang is currently offline  amin farhang
Messages: 39
Registered: November 2010
Member
Dear All,

I have a 16x16 sparse matrix which its values are big. Unfortunately IDL return wrong value for matrix inverse. Therefore when I run command (print,invert(A)##A) the returned is not an identity matrix. I check the singularity of matrix and this inversion retuned correctly by other softwares like MATLAB or MATHEMATICA or even FORTRAN.
What is happening?


EXAMPLE:

A = [ [7.3339770e12, 0.0, 0.0, 0.0, 0.0, 7.3339770e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] , $
[0.0, 5.4254596e12, 0.0, 0.0, 0.0, 5.4254596e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
[0.0, 0.0, 4.9832916e13, 0.0, 0.0, 0.0, 4.9832916e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0, 9.7295220e13, 0.0, 0.0, 9.7295220e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0, 0.0, 2.2853478e12, 2.2853478e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
[7.3339770e12, 5.4254596e12, 0.0, 0.0, 2.2853478e12, 1.5044784e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
[0.0, 0.0, 4.9832916e13, 9.7295220e13, 0.0, 0.0, 1.5037721e14, 3.2490665e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.2490665e12, 3.2490665e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7204194e13, 8.7204194e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7204194e13, 1.9651177e14, 0.0, 0.0, 1.3385090e13, 9.5922483e13, 0.0, 0.0], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.6135226e12, 1.1430911e12, 0.0, 0.0, 1.9939955e12, 1.4764362e12], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1430911e12, 1.1430911e12, 0.0, 0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.3385090e13, 0.0, 0.0, 1.3385090e13, 0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.5922483e13, 0.0, 0.0, 0.0, 9.5922483e13, 0.0, 0.0], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.9939955e12, 0.0, 0.0, 0.0, 1.9939955e12, 0.0], $
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.4764362e12, 0.0, 0.0, 0.0, 0.0, 1.4764362e12] ]


BUT IDL return wrong inverse matrix

IDL> B = invert(A,/double)
IDL> print, B

1.0e-05 *

[-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
[-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
[0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
[-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
[0.134, 0.134, 0.0, 0.0, 0.134,-0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
[0.0, 0.0,-0.0152,-0.0152, 0.0, 0.0, 0.0152,-0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
[ 0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0167, 0.0167, 0.0, 0.0,-0.0167,-0.0167, 0.0, 0.0], $
[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-1.342, 1.342,-0.0,-0.0, 1.342, 1.342], $
[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342], $
[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
[-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342], $
[-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342] ]


best regards,

Report message to a moderator

Re: IDL inverse matrix problem [message #89977 is a reply to message #89975] Mon, 12 January 2015 11:35 Go to previous messageGo to next message
Lajos Foldy is currently offline  Lajos Foldy
Messages: 176
Registered: December 2011
Senior Member
On Monday, January 12, 2015 at 6:56:36 PM UTC+1, Amin Farhang wrote:
> Dear All,
>
> I have a 16x16 sparse matrix which its values are big. Unfortunately IDL return wrong value for matrix inverse. Therefore when I run command (print,invert(A)##A) the returned is not an identity matrix. I check the singularity of matrix and this inversion retuned correctly by other softwares like MATLAB or MATHEMATICA or even FORTRAN.
> What is happening?
>
>
> EXAMPLE:
>
> A = [ [7.3339770e12, 0.0, 0.0, 0.0, 0.0, 7.3339770e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] , $
> [0.0, 5.4254596e12, 0.0, 0.0, 0.0, 5.4254596e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 4.9832916e13, 0.0, 0.0, 0.0, 4.9832916e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 9.7295220e13, 0.0, 0.0, 9.7295220e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 2.2853478e12, 2.2853478e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [7.3339770e12, 5.4254596e12, 0.0, 0.0, 2.2853478e12, 1.5044784e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 4.9832916e13, 9.7295220e13, 0.0, 0.0, 1.5037721e14, 3.2490665e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.2490665e12, 3.2490665e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7204194e13, 8.7204194e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7204194e13, 1.9651177e14, 0.0, 0.0, 1.3385090e13, 9.5922483e13, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.6135226e12, 1.1430911e12, 0.0, 0.0, 1.9939955e12, 1.4764362e12], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1430911e12, 1.1430911e12, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.3385090e13, 0.0, 0.0, 1.3385090e13, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.5922483e13, 0.0, 0.0, 0.0, 9.5922483e13, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.9939955e12, 0.0, 0.0, 0.0, 1.9939955e12, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.4764362e12, 0.0, 0.0, 0.0, 0.0, 1.4764362e12] ]
>
>
> BUT IDL return wrong inverse matrix
>
> IDL> B = invert(A,/double)
> IDL> print, B
>
> 1.0e-05 *
>
> [-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [0.134, 0.134, 0.0, 0.0, 0.134,-0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [0.0, 0.0,-0.0152,-0.0152, 0.0, 0.0, 0.0152,-0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0167, 0.0167, 0.0, 0.0,-0.0167,-0.0167, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-1.342, 1.342,-0.0,-0.0, 1.342, 1.342], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
> [-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342], $
> [-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342] ]
>
>
> best regards,

Scale your matrix, use double precision and LAPACK, eg:

m=double(max(abs(a)))
inverse=la_invert(a/m)/m

regards,
Lajos

Report message to a moderator

Re: IDL inverse matrix problem [message #89978 is a reply to message #89975] Mon, 12 January 2015 11:37 Go to previous message
wlandsman is currently offline  wlandsman
Messages: 743
Registered: June 2000
Senior Member
You want to specify your matrix as double precision and then you will get the expected answer. It is *not* enough to use the /DOUBLE keyword to INVERT(). See, for example,

http://www.idlcoyote.com/math_tips/sky_is_falling.html

--Wayne


On Monday, January 12, 2015 at 12:56:36 PM UTC-5, Amin Farhang wrote:
> Dear All,
>
> I have a 16x16 sparse matrix which its values are big. Unfortunately IDL return wrong value for matrix inverse. Therefore when I run command (print,invert(A)##A) the returned is not an identity matrix. I check the singularity of matrix and this inversion retuned correctly by other softwares like MATLAB or MATHEMATICA or even FORTRAN.
> What is happening?
>
>
> EXAMPLE:
>
> A = [ [7.3339770e12, 0.0, 0.0, 0.0, 0.0, 7.3339770e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] , $
> [0.0, 5.4254596e12, 0.0, 0.0, 0.0, 5.4254596e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 4.9832916e13, 0.0, 0.0, 0.0, 4.9832916e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 9.7295220e13, 0.0, 0.0, 9.7295220e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 2.2853478e12, 2.2853478e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [7.3339770e12, 5.4254596e12, 0.0, 0.0, 2.2853478e12, 1.5044784e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 4.9832916e13, 9.7295220e13, 0.0, 0.0, 1.5037721e14, 3.2490665e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.2490665e12, 3.2490665e12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7204194e13, 8.7204194e13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7204194e13, 1.9651177e14, 0.0, 0.0, 1.3385090e13, 9.5922483e13, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.6135226e12, 1.1430911e12, 0.0, 0.0, 1.9939955e12, 1.4764362e12], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1430911e12, 1.1430911e12, 0.0, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.3385090e13, 0.0, 0.0, 1.3385090e13, 0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.5922483e13, 0.0, 0.0, 0.0, 9.5922483e13, 0.0, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.9939955e12, 0.0, 0.0, 0.0, 1.9939955e12, 0.0], $
> [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.4764362e12, 0.0, 0.0, 0.0, 0.0, 1.4764362e12] ]
>
>
> BUT IDL return wrong inverse matrix
>
> IDL> B = invert(A,/double)
> IDL> print, B
>
> 1.0e-05 *
>
> [-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [-0.134,-0.134, 0.0, 0.0,-0.134, 0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [0.134, 0.134, 0.0, 0.0, 0.134,-0.134, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [0.0, 0.0,-0.0152,-0.0152, 0.0, 0.0, 0.0152,-0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0152, 0.0152, 0.0, 0.0,-0.0152, 0.0152, 0.0, 0.0, 0.0, 0.0,-0.0,-0.0, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.0167, 0.0167, 0.0, 0.0,-0.0167,-0.0167, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-1.342, 1.342,-0.0,-0.0, 1.342, 1.342], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
> [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0167,-0.0167, 0.0, 0.0, 0.0167, 0.0167, 0.0, 0.0], $
> [-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342], $
> [-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0, 1.342,-1.342, 0.0, 0.0,-1.342,-1.342] ]
>
>
> best regards,

Report message to a moderator

  Switch to threaded view of this topic Create a new topic Submit Reply
Previous Topic: Re: FG subplot question
Next Topic: For loop - Question

-=] Back to Top [=-
[ Syndicate this forum (XML) ] [ RSS ] [ PDF ]

Current Time: Wed Oct 08 13:40:03 PDT 2025

Total time taken to generate the page: 0.00470 seconds
.:: Contact :: Home ::.

Powered by: FUDforum 3.0.7.
Copyright ©2001-2013 FUDforum Bulletin Board Software