[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