;$Id: la_determ.pro,v 1.5 2004/01/21 15:54:55 scottm Exp $
;
; Copyright (c) 2001-2004, Research Systems, Inc.  All rights reserved.
;   Unauthorized reproduction prohibited.
;+
; NAME:
;   LA_DETERM
;
; PURPOSE:
;   This function computes the determinant of an N by N array.
;
; CALLING SEQUENCE:
;   Result = LA_DETERM(A)
;
; INPUTS:
;   A: An N by N array of any numeric type. A may be complex.
;
; KEYWORD PARAMETERS:
;   CHECK:  If set to a non-zero value, A is checked for singularity.
;           The determinant of a singular array is returned as zero if
;           this keyword is set. Run-time errors may result if A is
;           singular and this keyword is not set.
;
;   DOUBLE: If set to a non-zero value, computations are done in
;           double precision arithmetic.
;
;   ZERO:   Use this keyword to set the value of floating-point
;           zero. A floating-point zero on the main diagonal of
;           a triangular matrix results in a zero determinant.
;           For single-precision inputs, the default value is
;           1.0e-6. For double-precision inputs, the default value
;           is 1.0e-12.
;
; EXAMPLE:
;       Define an array (a).
;         a = [[ 2.0,  1.0,  1.0], $
;              [ 4.0, -6.0,  0.0], $
;              [-2.0,  7.0,  2.0]]
;       Compute the determinant.
;         result = LA_DETERM(a)
;
; PROCEDURE:
;       LU decomposition is used to represent the input array in
;       triangular form. The determinant is computed as the product
;       of diagonal elements of the triangular form. Row interchanges
;       are tracked during the LU decomposition to ensure the correct
;       sign, + or - .
;
; REFERENCE:
;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
;       Erwin Kreyszig
;       ISBN 0-471-55380-8
;
;       Anderson et al., LAPACK Users' Guide, 3rd ed., SIAM, 1999.
;
; MODIFICATION HISTORY:
;       Written by:  CT, RSI, December 2001. Similar to determ.pro but
;            uses LAPACK LA_LUDC for complex input.
;-

FUNCTION la_determ, A, $
    Check = Check, $
    Double = DoubleIn, $
    Zero = ZeroIn

    compile_opt idl2

    ON_ERROR, 2  ;Return to caller if error occurs.

    dims = SIZE(A, /DIMENSION)
    if (N_ELEMENTS(dims) ne 2) then dims = [0,-1]
    if (dims[0] ne dims[1]) then $
        MESSAGE, 'Input must be a square array.'

    type = SIZE(A, /TYPE)
    isComplex = (type eq 6) or (type eq 9)

    ; If the DOUBLE keyword is not set then the internal precision and
    ; result are identical to the type of input.
    Double = (N_ELEMENTS(DoubleIn) gt 0) ? KEYWORD_SET(DoubleIn) : $
        ((type eq 5) or (type eq 9))
    zeroRet = isComplex ? COMPLEX(0, DOUBLE=Double) : (Double ? 0d : 0.0)


    ; If ZERO keyword is not set then use the default zero cutoff.
    Zero = (N_ELEMENTS(ZeroIn) gt 0) ? $
        ABS(ZeroIn[0]) : (Double ? 1d-12 : 1e-6)

    ; Make a copy of the array for its LU decomposition.
    ALUd = A

    ; Compute LU decomposition, either with or without STATUS.
    if KEYWORD_SET(Check) then begin
        LA_LUDC, ALUd, Index, $
            Double = Double, $
            Interchanges = Sign, $
            STATUS=status
        ; Singular matrix? Return zero.
        if (status ne 0) then $
            RETURN, zeroRet
    endif else begin
        ; This will throw an error for a singular matrix.
        LA_LUDC, ALUd, Index, $
            Double = Double, $
            Interchanges = Sign
    endelse

    ; We only need the main diagonal.
    ALUd = DIAG_MATRIX(ALUd)

    ;Are there any zeros on the main diagonal?
    Cnt = TOTAL( ABS(ALUd) le Zero)
    if Cnt ne 0 then $ ;A zero on the main diagonal results in a zero determ.
        RETURN, zeroRet

    Det = PRODUCT(TEMPORARY(ALUd), DOUBLE=Double)

    RETURN, Sign * Det

end