
determinant
Calculate the determinant of a small square real matrix using a partial-pivoting Gaussian elimination scheme.
Available in version 6.3.0 and later.
Prototype
function determinant ( x [*][*] : numeric ) return_val [1] : float or double
Arguments
xA two-dimensional variable of numeric type. The dimension sizes must be the same (i.e. [N][N]). N should be 'small' (say, .le. 7).
Return value
A scalar of type double if x is type double; otherwise the returned value will be type float.
Description
Normally, NCL would use an LAPACK code to compute a mathematical quantity such as the determinant. However, LAPACK contains no such subroutine due to issues with accuracy and stability:
http://www.netlib.org/lapack/faq.html#_are_there_routines_in_lapack_to_compute_determinantsBasically, numerical analysts do not want people calculating the determinant. :-(
See Also
WARNING: User Beware: This subroutine has been tested using small array sizes only.
See examples.
An Introduction to Computational Physics
Tao Pang
Cambridge University Press, 1997
Specifically, the determinant of the square (NxN) matrix is calculated using
the partial-pivoting Gaussian elimination scheme. A more robust approach would
be to use (say) QR decomposition.
Examples
The following use integer argument for convenience. All values sent to the determinant function are promoted to double prior to the computation. The computations are done in double precision.
Example 1:
a = (/ (/ 1, 0, 2, -1/) \ ; N=4 , (/ 3, 0, 0, 5/) \ , (/ 2, 1, 4, -3/) \ , (/ 1, 0, 5, 0/) /) d4 = determinant(a) ; d4=30.0Example 3: The following matches the result from the online determinant calculator.
b = (/ (/ 1,1,0,0,0,0 /) \ ; N=6 , (/ 2,2,2,0,0,0 /) \ , (/ 0,3,3,3,0,0 /) \ , (/ 0,0,4,4,4,0 /) \ , (/ 0,0,0,5,5,5 /) \ , (/ 0,0,0,0,6,6 /) /) d6 = determinant(b) ; d6=724.0
Example 4: The following matches the result from the online determinant calculator.
c = (/ (/ 1,1,0,0,0,0,0, /) \ ; N=7 , (/ 2,2,2,0,0,0,0, /) \ , (/ 0,3,3,3,0,0,0, /) \ , (/ 0,0,4,4,4,0,0, /) \ , (/ 0,0,0,5,5,5,0, /) \ , (/ 0,0,0,0,6,6,6, /) \ , (/ 0,0,0,0,0,7,7, /) /) d7 = determinant(c) ; d7=5040
Example 5: The following matches the result from the online determinant calculator.
e = (/ (/ 1,1,0,0,0,0,0,0, /) \ ; N=8 , (/ 2,2,2,0,0,0,0,0, /) \ , (/ 0,3,3,3,0,0,0,0, /) \ , (/ 0,0,4,4,4,0,0,0, /) \ , (/ 0,0,0,5,5,5,0,0, /) \ , (/ 0,0,0,0,6,6,6,0, /) \ , (/ 0,0,0,0,0,7,7,7, /) \ , (/ 0,0,0,0,0,0,8,8, /) /) d8 = determinant(e) ; d8=0
Example 6:
f = (/ (/ 1, 0.5, 0, 0, 0, 0, 0, 0, 0 /) \ ; N=9 , (/-1, 0.5, 1, 0, 0, 0, 0, 0, 0 /) \ , (/ 0,-1, 0, 2.5, 0, 0, 0, 0, 0 /) \ , (/ 0, 0, -1,-1.5, 4, 0, 0, 0, 0 /) \ , (/ 0, 0, 0,-1, -3, 5.5, 0, 0, 0 /) \ , (/ 0, 0, 0, 0, -1,-4.5, 7, 0, 0 /) \ , (/ 0, 0, 0, 0, 0,-1, -6, 8.5, 0 /) \ , (/ 0, 0, 0, 0, 0, 0, -1,-7.5,10 /) \ , (/ 0, 0, 0, 0, 0, 0, 0,-1, -9 /) /) d9 = determinant(f) ; d9=1.0 This is correct.
Example 7: The following matches the result from the online determinant calculator.
e = (/ (/ 1,1,0,0,0,0,0,0,0,0 /) \ ; N=10 , (/ 2,2,2,0,0,0,0,0,0,0 /) \ , (/ 0,3,3,3,0,0,0,0,0,0 /) \ , (/ 0,0,4,4,4,0,0,0,0,0 /) \ , (/ 0,0,0,5,5,5,0,0,0,0 /) \ , (/ 0,0,0,0,6,6,6,0,0,0 /) \ , (/ 0,0,0,0,0,7,7,7,0,0 /) \ , (/ 0,0,0,0,0,0,8,8,8,0 /) \ , (/ 0,0,0,0,0,0,0,9,9,9 /) \ , (/ 0,0,0,0,0,0,0,0,10,10 /) /) d10 = determinant(e) ; d10=-3628800.0 ????????