NCL Home >
Documentation >
Functions >
General applied math
inverse_matrix
Computes the inverse of a square matrix using LU factorization.
Prototype
function inverse_matrix ( A [*][*] : numeric ) return_val [dimsizes(A)] : float or double
Arguments
AA two-dimensional (N, N) array.
Return value
An array of the same size, shape and type A.
Description
NOTE: NCL does not have a numerical analysis focus. There are better tools for inverting matrices: Matlab, Octave, IDL and Python.
NOTE: Generally, it is best not to invert a matrix to solve an equation. See: Do not invert that matrix.
Perhaps, NCL's solve_linsys would be appropriate.
Using the LAPACK routines DGETRF and DGETRI, this function computes the inverse of a general matrix A(N, M) [here: N=M]using the LU factorization method. This method inverts U and then computes inv(A) by solving the system inv(A) * L = inv(U) for inv(A).
The returned inverse provides one way of solving the set of equations A * x = b (i.e. x = inv(A) * b). Thus, if the inverse is known, the vector x can be solved for any b. However, if efficiency is a concern, this is not the preferred method (see solve_linsys).
See Also
Examples
The identity matrix (I) is found via I = aInv # a. The "#" operator computes the dot product of two arrays.
Example 1: Compute the determinant and inverse of a 3 x 3 matrix:
a = (/(/1.0, -1.0, 2.0/) \
,(/3.0, 0.0, 1.0/) \
,(/1.0, 0.0, 2.0/) /)
write_matrix(a, 3f4.1, False)
aDet = determinant(a)
aInv = inverse_matrix(a) ; [3] x [3]
aId = aInv#a
print("aDet="+aDet)
write_matrix(aInv, "3f4.1", False)
write_matrix(aId , "3f4.1", False)
----------------
qDet=5
1.0 -1.0 2.0 <=== original matrix
3.0 0.0 1.0
1.0 0.0 2.0
===
0.0 0.4 -0.2 <=== inverse
-1.0 -0.0 1.0
0.0 -0.2 0.6
===
1.0 0.0 0.0 <=== identity
0.0 1.0 0.0
0.0 0.0 1.0
Example 2: Compute the determinant and inverse of a 4 x 4 matrix:
x = (/ (/3.0, 2.0, -1.0, 2.0/), \
(/1.0, 4.0, 0.0, 2.0/), \
(/2.0, 1.0, 2.0, -1.0/), \
(/1.0, 1.0, -1.0, 3.0/) /)
write_matrix(x, 4f4.1, False)
xDet = determinant(x)
xInv = inverse_matrix(x) ; [4] x [4]
xId = xInv#x
print("xDet="+xDet)
write_matrix(xInv, "4f4.1", False)
write_matrix(toint(xId) , "4i4" , False)
----------------
xDet=44
3.0 2.0 -1.0 2.0
1.0 4.0 0.0 2.0
2.0 1.0 2.0 -1.0
1.0 1.0 -1.0 3.0
===
0.3 -0.2 0.1 -0.0
0.1 0.3 -0.1 -0.3
-0.5 0.0 0.5 0.5
-0.3 -0.0 0.2 0.6
===
1.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0
0.0 0.0 1.0 0.0
0.0 0.0 0.0 1.0
Example 3: Compute the inverse of a NxN matrix (N=20). Due to computer numerics,
values such as -0.0 and 0.99999 are returned. These are manipulated
using round and toint to return the expected value [0.0 and 1.0].
b = random_normal( 0,10, (/N,N/)) write_matrix(b, "20f6.1", False) bDet = determinant(b) bInv = inverse_matrix(b) bId = bInv#b print("bDet="+bDet) bId = where(abs(bId).lt.1e-5, 0.0, bId) ; cleaner output write_matrix(toint(round(bID,0),"20i3",False) ; printing convenience ----------------- bDet=-6.023e+28 8.7 8.7 19.5 9.8 3.5 -1.3 -8.6 8.8 9.3 -11.7 9.2 7.4 -28.2 0.1 0.5 7.5 1.8 -6.2 10.2 -5.0 11.7 -6.7 -3.3 16.9 15.1 0.3 3.5 1.7 -13.8 8.1 20.8 9.1 -11.3 -2.9 17.3 3.2 14.2 18.0 -5.2 1.6 -6.2 -0.5 -4.1 3.5 13.9 -0.8 -12.3 11.1 3.1 -5.4 -28.3 -11.8 15.9 -1.9 -11.1 20.9 -16.3 -11.3 5.6 -6.6 4.2 -6.1 6.0 -5.0 11.1 13.9 -1.2 10.6 18.6 2.2 3.1 1.3 3.2 15.5 -9.8 -13.7 -1.6 1.8 10.8 -11.3 9.9 -7.2 -1.1 10.3 8.0 7.1 -2.6 -7.7 5.7 8.7 4.3 16.8 9.5 -6.0 -3.5 1.7 -5.0 11.6 -4.6 -0.0 -9.0 -3.2 -2.0 -16.4 4.8 2.2 -3.0 -6.9 8.5 15.7 -0.9 7.6 29.0 5.1 -1.6 -10.0 6.8 -11.9 -3.7 -4.8 -4.1 6.2 -14.3 11.4 -10.9 0.9 -3.4 7.2 10.2 3.7 -0.1 4.3 12.5 -11.2 11.9 6.7 2.5 4.8 12.5 9.0 -5.6 -7.2 5.4 9.7 -22.7 2.7 5.6 7.6 4.8 -9.7 1.7 -6.2 24.3 -7.8 -4.9 -3.3 1.0 -21.3 -18.4 17.8 -10.5 16.5 19.0 -6.5 -14.1 -2.3 13.7 -0.4 8.5 -3.2 3.0 -23.0 7.0 -5.5 20.9 19.1 0.1 -12.8 -6.7 -5.5 9.8 8.8 17.0 -1.0 -4.1 3.7 17.7 3.6 22.6 -25.5 -13.8 8.1 -2.8 11.1 -3.0 0.1 7.9 5.4 9.3 15.2 -3.9 4.1 -10.6 10.8 3.4 7.3 -13.2 6.7 -11.0 11.6 3.0 14.7 0.2 -5.8 -22.6 -8.8 -18.7 10.9 -6.1 0.5 -7.4 11.4 -3.8 -21.2 -6.0 -9.6 -4.3 5.4 9.4 5.4 -7.1 4.2 6.8 -7.8 12.3 -0.8 -3.5 -1.6 -11.6 -3.8 2.0 16.6 -2.9 -6.6 10.9 5.7 2.8 12.2 4.9 4.6 11.2 -10.3 -9.7 -19.4 1.3 -14.0 -14.1 8.9 8.1 7.1 4.7 1.4 13.6 -8.6 -12.7 -1.1 -1.9 -1.8 7.4 -7.2 -4.0 -7.4 3.8 -6.3 -4.9 6.2 -5.2 -14.1 10.0 -0.5 9.8 -13.3 -3.5 9.7 -6.9 -10.0 -11.4 6.2 -3.9 1.0 -10.4 -7.5 17.6 1.9 -8.9 -0.3 -1.1 13.4 -8.4 6.0 12.5 -1.7 -0.3 2.0 -0.8 -2.2 4.4 -5.1 -2.6 -14.9 0.2 13.9 5.7 2.4 -13.7 -0.2 -9.4 0.0 7.1 6.5 4.4 11.5 -10.8 10.7 5.6 12.0 -9.2 2.7 9.2 8.8 0.7 9.5 7.6 3.2 5.8 -6.1 -7.9 7.1 2.3 -5.0 -1.9 -7.9 4.8 7.9 -19.1 7.7 3.0 -11.7 -7.8 -7.5 -6.1 -7.0 -15.2 24.6 -0.8 -8.1 -17.4 -11.2 -4.6 5.0 4.3 -5.7 7.2 3.0 2.9 -3.0 3.6 -15.1 -14.3 6.6 10.8 -3.5 -12.7 -3.2 -3.7 -14.5 0.3 -8.8 16.2 -19.0 -3.1 -2.4 -1.7 -8.2 -11.0 24.6 -9.3 -3.8 9.3 -5.2 -6.6 -9.9 19.2 -0.2 16.6 15.4 -7.3 -16.3 -3.9 -7.9 -0.0450 0.0339 0.0001 0.0072 -0.0108 -0.0399 -0.0412 0.0143 -0.0226 -0.0149 -0.0232 0.0040 -0.0220 0.1120 -0.0261 -0.0096 0.0795 -0.0107 -0.0256 -0.0314 -0.0055 -0.0049 -0.0075 -0.0195 -0.0247 0.0011 -0.0207 -0.0025 0.0088 0.0031 0.0059 -0.0152 -0.0048 0.0213 -0.0061 0.0000 0.0503 -0.0139 -0.0107 -0.0106 0.0304 0.0221 0.0084 -0.0108 0.0034 0.0467 0.0114 -0.0118 0.0194 0.0372 0.0350 0.0072 0.0187 -0.0222 0.0322 -0.0067 -0.0483 0.0176 0.0275 0.0206 0.0106 -0.0091 0.0070 -0.0059 0.0111 -0.0037 0.0075 0.0114 0.0063 0.0027 -0.0041 -0.0093 -0.0024 -0.0255 0.0058 -0.0058 0.0120 0.0020 0.0143 -0.0059 0.0190 -0.0032 0.0147 -0.0109 0.0089 0.0243 -0.0068 -0.0031 0.0003 0.0018 -0.0129 -0.0097 0.0220 -0.0397 0.0055 0.0088 -0.0116 0.0013 0.0019 0.0049 -0.0165 0.0147 -0.0030 0.0002 -0.0076 -0.0061 -0.0100 0.0015 -0.0095 0.0175 0.0214 -0.0112 -0.0027 0.0241 -0.0141 -0.0219 0.0064 -0.0052 -0.0008 0.0205 -0.0522 0.0244 0.0104 0.0258 0.0046 -0.0398 -0.0092 0.0164 0.0027 -0.0002 -0.0038 0.0132 -0.0223 0.0504 -0.0434 -0.0056 0.0241 -0.0025 -0.0047 -0.0320 -0.0199 0.0482 0.0157 0.0310 -0.0304 -0.0156 0.0033 0.0137 -0.0049 0.0057 0.0215 0.0273 -0.0129 0.0457 -0.0139 -0.0065 0.0043 0.0013 -0.0039 -0.0112 0.0083 -0.0431 -0.0094 0.0127 0.0373 -0.0193 0.0013 0.0064 -0.0020 -0.0225 -0.0369 -0.0064 -0.0001 -0.0144 -0.0077 -0.0051 0.0122 -0.0057 -0.0179 -0.0152 -0.0434 0.0397 0.0089 0.0107 0.0027 -0.0100 -0.0108 0.0024 -0.0094 -0.0096 0.0045 0.0014 -0.0209 0.1093 -0.0348 -0.0420 0.0420 0.0104 -0.0262 -0.0436 0.0223 -0.0198 -0.0207 0.0109 0.0019 0.0086 0.0077 -0.0036 0.0183 -0.0308 -0.0094 -0.0194 0.0100 -0.0328 0.0073 0.0353 -0.0137 -0.0050 -0.0153 0.0037 0.0035 0.0232 0.0067 0.0023 0.0015 0.0096 0.0127 0.0008 -0.0050 0.0158 0.0218 0.0260 -0.0109 0.0031 -0.0113 0.0074 -0.0150 0.0069 0.0103 0.0030 0.0079 0.0046 0.0024 -0.0027 -0.0083 0.0274 0.0100 -0.0030 0.0165 0.0096 0.0106 -0.0024 0.0082 -0.0268 0.0234 0.0222 -0.0120 -0.0018 0.0132 0.0143 -0.0013 0.0079 -0.0025 0.0108 -0.0189 0.0117 -0.0032 -0.0102 0.0090 -0.0097 -0.0034 -0.0106 -0.0095 0.0090 0.0008 0.0128 0.0156 0.0135 -0.0217 -0.0120 0.0181 0.0237 0.0108 -0.0073 -0.0014 0.0245 0.0211 -0.0053 0.0036 0.0172 -0.0001 0.0251 0.0182 -0.0342 0.0178 0.0065 -0.0253 0.0214 0.0232 0.0213 -0.0165 0.0085 0.0019 0.0077 0.0047 -0.0180 -0.0042 -0.0075 0.0061 -0.0196 -0.0018 -0.0127 -0.0209 0.0459 -0.0237 0.0009 0.0101 -0.0060 -0.0339 -0.0148 -0.0081 -0.0041 -0.0122 -0.0093 -0.0178 -0.0066 -0.0139 0.0055 -0.0135 0.0063 -0.0036 -0.0185 -0.0178 0.0189 -0.0060 -0.0247 0.0229 -0.0255 -0.0051 -0.0018 0.0152 -0.0172 -0.0128 -0.0004 0.0101 0.0094 0.0177 -0.0231 0.0188 0.0167 0.0190 -0.0018 0.0131 -0.0566 0.0360 0.0077 -0.0411 -0.0016 0.0133 0.0237 0.0055 0.0116 0.0063 0.0030 -0.0120 0.0212 0.0313 -0.0158 0.0022 0.0109 0.0109 -0.0050 0.0087 0.0026 0.0068 0.0082 -0.0237 0.0093 0.0138 0.0061 0.0083 0.0126 0.0024 -0.0211 -0.0006 0.0338 0.0072 -0.0056 -0.0090 0.0048 -0.0024 -0.0182 0.0231 0.0161 0.0028 -0.0047 -0.0127 0.0197 -0.0243 -0.0050 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Errors
The LAPACK subrourines return the following information:
- info
- = 0: successful exit
- < 0: if info = -i, the i-th argument had an illegal value
- > 0: if info = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.
- < 0: if info = -i, the i-th argument had an illegal value