[ncl-talk] inverse_matrix problem
Walter Hannah
whannah at atmos.colostate.edu
Mon Nov 10 16:25:15 MST 2014
Will,
I tried your suggestion, but all I got back was missing values.
After Dennis's suggestion, I tried replacing the zeros with tiny random
numbers, using this code:
R = random_normal( 0.,1e-12, (/nlev*2,nlev*2/))
M = where(M.eq.0,R,M)
iM = inverse_matrix(M)
Interestingly, this worked perfectly, and I got what I expected.
I also tried using this approach to trying Will's suggestion, and that
also worked with the random numbers!
I can't figure out a mathematical explanation, but this seems like a
pretty nice fix. Might be good to add this to the online documentation.
Thanks,
Walter
On 11/10/14 6:17 PM, Will Hobbs wrote:
> Walter
>
> Instead of using the inverse matrix() function, try using
> solve_linsys(A, B) instead, for A = your matrix and B = I (the
> identity matrix). Since A * inv(A) = I, the output will be the inverse.
>
> Mathematically this is no different to inverse_matrix(), but the code
> is a little different. In particular, if the matrix is singular (i.e.
> the determinant = 0, and it cannot be inverted) then it will return an
> error code.
> https://www.ncl.ucar.edu/Document/Functions/Built-in/solve_linsys.shtml
>
> This might help you figure out if the matrix itself is the problem (I
> suspect, as Dennis has suggested, that may be the issue).
>
> Will
>
>
> From: Walter Hannah <whannah at atmos.colostate.edu
> <mailto:whannah at atmos.colostate.edu>>
> Date: Tuesday, 11 November 2014 9:52 AM
> To: Dennis Shea <shea at ucar.edu <mailto:shea at ucar.edu>>
> Cc: "ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>" <ncl-talk at ucar.edu
> <mailto:ncl-talk at ucar.edu>>
> Subject: Re: [ncl-talk] inverse_matrix problem
>
> Dennis,
>
> Wow! Thanks for going above and beyond! This helps a lot. I'll try to
> look for a way to get a modified matrix that doesn't have all those zeros.
>
> Walter
>
> On 11/10/14 5:50 PM, Dennis Shea wrote:
>> You have rows of 0.0 ... I think that is causing issues. I'm sure
>> some 'math-whiz' could answer but I attached a simple test code that
>> reads your data (thx for providing). It calculates the inverse and
>> the identity matrix. Note that NCL is using a direct call to the
>> documented LAPACK routines.
>>
>> ===
>> TEST1: duplicates your findings
>>
>> TEST 2: create a bogus 20x20 array of random numbers. No row consists
>> of all constants. Identity matrix is successfully calculated.
>>
>> TEST 3: wherever the original array had a row of all 0.0, substitute
>> random numbers. The identity matrix is successfully calculated.
>>
>> On Mon, Nov 10, 2014 at 1:57 PM, Walter Hannah
>> <whannah at atmos.colostate.edu <mailto:whannah at atmos.colostate.edu>> wrote:
>>
>> Hi,
>>
>> I'm trying to use the "inverse_matrix" function. The matrix
>> doesn't have any missing values, but I still get this error message:
>>
>> Variable: M
>> Type: float
>> Total Size: 1600 bytes
>> 400 values
>> Number of Dimensions: 2
>> Dimensions and sizes: [pertlev | 20] x [tendlev | 20]
>> Coordinates:
>> pertlev: [50..1950]
>> tendlev: [50..1950]
>> Number Of Attributes: 2
>> missing_value : 9.96921e+36
>> _FillValue : 9.96921e+36
>> (0)
>> ----------------------------------------------------------------------------
>>
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00
>> 0.00 -0.31 -0.06 0.15 -0.27 -0.56 -1.21 -1.38
>> -1.36 1.49 0.00 0.00 4.65 0.37 0.46 0.47 0.52
>> 0.55 0.77 2.49
>> 0.00 0.26 -0.35 -0.29 -0.65 -1.13 -2.24 -2.45
>> -2.47 2.91 0.00 0.00 4.07 2.45 0.95 0.87 0.90
>> 1.05 1.55 4.81
>> 0.00 -0.11 1.43 -0.91 -1.69 -1.52 -2.63 -2.81
>> -2.69 3.17 0.00 0.00 -2.35 2.61 1.49 1.01
>> 1.00 1.18 1.78 5.32
>> 0.00 -0.02 0.55 1.08 -1.74 -2.04 -2.59 -2.71
>> -2.39 2.63 0.00 0.00 -0.75 -1.10 1.35 1.02
>> 0.89 1.00 1.57 4.87
>> 0.00 -0.06 0.16 0.52 0.72 -1.86 -2.88 -2.74
>> -2.29 2.31 0.00 0.00 0.31 -1.30 -0.26 0.87 0.75
>> 0.82 1.27 4.50
>> 0.00 -0.13 -0.07 0.20 0.54 0.27 -2.74 -3.10
>> -2.27 2.04 0.00 0.00 0.08 -0.49 -0.47 0.14 0.72
>> 0.69 1.13 4.24
>> 0.00 -0.15 -0.06 0.20 -0.17 0.02 -0.77 -3.87
>> -3.13 2.39 0.00 0.00 0.69 -0.08 0.09 0.03 0.36
>> 0.90 1.29 4.59
>> 0.00 -0.14 0.05 0.09 -0.17 -0.29 -0.14 -0.66
>> -5.36 2.00 0.00 0.00 0.39 0.17 0.14 0.07 -0.01
>> 0.22 1.00 4.05
>> 0.00 0.02 -0.07 -0.01 0.00 0.02 0.43 0.94
>> 0.77 -2.90 0.00 0.00 -0.09 0.09 -0.04 -0.07 -0.17
>> -0.26 -0.49 0.54
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00
>> 0.00 -0.01 -0.04 0.12 -0.02 -0.01 -0.03 -0.05
>> -0.05 0.04 0.00 0.00 -2.17 -0.25 -0.14 0.01
>> 0.02 0.02 0.02 0.12
>> 0.00 0.00 -0.19 0.11 0.28 0.06 0.11 0.06
>> 0.07 -0.13 0.00 0.00 2.75 -1.97 -0.04 -0.12 -0.05
>> -0.05 -0.08 -0.10
>> 0.00 -0.02 -0.14 -0.08 0.31 0.21 0.13 0.08
>> 0.06 -0.08 0.00 0.00 0.48 1.11 -1.09 -0.13 -0.10
>> -0.06 -0.09 -0.12
>> 0.00 -0.01 -0.03 -0.11 -0.07 0.32 0.27 0.12
>> 0.18 -0.17 0.00 0.00 0.02 0.63 0.52 -0.94 -0.02
>> -0.15 -0.08 -0.25
>> 0.00 0.02 0.08 -0.06 -0.24 0.10 0.43 0.26
>> 0.16 -0.17 0.00 0.00 0.07 0.31 0.35 0.40 -1.06
>> -0.06 -0.20 -0.26
>> 0.00 0.06 0.03 -0.10 0.07 0.11 1.03 1.38
>> 0.53 -1.10 0.00 0.00 -0.44 0.14 -0.03 0.08 0.19
>> -1.61 -0.76 -1.16
>> 0.00 0.11 -0.03 -0.08 0.07 0.16 0.46 1.28
>> 1.65 -1.21 0.00 0.00 -0.25 0.07 -0.06 0.00 0.14
>> 0.41 -3.04 -0.66
>> 0.00 0.12 -0.30 -0.15 0.80 1.58 2.86 3.62
>> 5.02 -3.15 0.00 0.00 -1.23 -0.91 -0.80 -0.94 -0.85
>> -0.66 0.21 -10.64
>>
>> (0)
>> ----------------------------------------------------------------------------
>> (0) # missing values: 0
>> (0) # values = 0.0 : 111
>> warning:inverse_matrix: info = 1; missing values not allowed
>>
>>
>> I thought it might have something to do with those zero values,
>> but if I replace them with non-zero values I get a similar, but
>> not identical, result:
>>
>> ...
>> (0)
>> ----------------------------------------------------------------------------
>> (0) # missing values: 0
>> (0) # values = 0.0 : 0
>> warning:inverse_matrix: info = 20; missing values not allowed
>>
>>
>>
>> Thanks,
>> Walter
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20141110/e8810619/attachment.html
More information about the ncl-talk
mailing list