# Source code for cond.pro:

```;\$Id: cond.pro,v 1.19 2001/05/23 21:20:32 chris Exp \$
;
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       COND
;
; PURPOSE:
;       This function computes the condition number of an N by N array.
;
; CATEGORY:
;       Complex Linear Algebra.
;
; CALLING SEQUENCE:
;       Result = COND(A)
;
; INPUTS:
;       A:      An N by N real or complex array.
;
; KEYWORD PARAMETERS:
;       DOUBLE: If set to a non-zero value, computations are done in
;               double precision arithmetic.
;
;   LNORM: Set this keyword to indicate which norm to use for the computation.
;          The possible values are:
;           LNORM=0  Use the L(Infinity) norm (maximum absolute row sum norm).
;           LNORM=1  Use the L(1) norm (maximum absolute column sum norm).
;                    For LNORM=0 or 1 the array A must be square.
;           LNORM=2  Use the L(2) norm (spectral norm), defined as
;                    the largest singular value, computed from SVDC.
;                    For LNORM=2 the array A cannot be complex.
;           If LNORM is not specified then LNORM=0 is used.
;
; EXAMPLE:
;       Define a complex array (a).
;         a = [[complex(1, 0), complex(2,-2), complex(-3,1)], \$
;              [complex(1,-2), complex(2, 2), complex(1, 0)], \$
;              [complex(1, 1), complex(0, 1), complex(1, 5)]]
;       Compute the condition number of the complex array (a) using
;       double-precision complex arithmetic.
;         result = cond(a, /double)
;
; PROCEDURE:
;    This function returns the condition number of an N x N real or
;    complex array A.
;       For LNORM=0 or 1, the condition number is norm(A)*norm(A_inverse).
;       If A is real and A_inverse is invalid (due to the singularity of A
;       or floating-point errors in the invert function), the condition
;       number is returned as a -1. If A is complex and A_inverse is invalid
;       (due to the singularity of A), calling COND results in floating-
;       point errors.
;       For LNORM=2, the condition number is defined as the ratio of the
;       largest to smallest singular values, computed using SVDC.
;       If the smallest singular value is zero, then Infinity is returned.
;       For LNORM=2 the array A cannot be complex.
;
; REFERENCE:
;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
;       Erwin Kreyszig
;       ISBN 0-471-55380-8
;
;       CRC Concise Encyclopedia of Mathematics
;       Eric W. Weisstein
;       ISBN 0-8493-9640-9
;
; MODIFICATION HISTORY:
;       Written by:  GGS, RSI, April 1992
;       Modified:    GGS, RSI, February 1994
;                    Accepts complex inputs. Added DOUBLE keyword.
;       Modified:    GGS, RSI, November 1994
;                    Added support for double-precision complex inputs.
;                    Changed NR_INVERT to INVERT.
;       Modified:    GGS, RSI, April 1996
;                    Modified keyword checking and use of double precision.
;       Modified: CT, RSI, May 2001
;             Allows L(1), L(2), L(Infinity) norm to be used. L(2) uses SVDC.
;-

function Cond, Array, \$
DOUBLE=double, \$
LNORM=lnormIn

COMPILE_OPT idl2

type = SIZE(Array, /TYPE)
dim = SIZE(Array, /DIMENSIONS)
if (N_ELEMENTS(dim) ne 2) then \$
MESSAGE, 'Input must be a two-dimensional array.'

dbl = (N_ELEMENTS(double) gt 0) ? KEYWORD_SET(double) : \$
(type eq 5) or (type eq 9)

; Default is to use L(Infinity) norm
lnorm = (N_ELEMENTS(lnormIn) gt 0) ? FIX(lnormIn[0]) : 0
if (lnorm lt 0) or (lnorm gt 2) then \$
MESSAGE, 'Keyword LNORM must be equal to 0, 1, or 2'

; For L(2) we need to use SVDC, otherwise call NORM.
if (lnorm eq 2) then begin

is_complex = (type eq 6) or (type eq 9)
if is_complex then \$
MESSAGE, 'Complex input not allowed for LNORM=2'

; If A is an M x N array, and M > N, then SVDC will return M
; singular values, with (M - N) of them set to garbage (zero).
; To avoid this, set the /COLUMN keyword if M > N.
SVDC, Array, W, U, V, \$
COLUMN=(dim[0] gt dim[1]), \$
DOUBLE=dbl
minSV = MIN(W, MAX=maxSV)

; Is the matrix singular? If so return Infinity without
; throwing an overflow message.
if (minSV eq 0) then \$
return, dbl ? !VALUES.d_infinity : !VALUES.f_infinity

; If the matrix is ill-conditioned this could still overflow.
return, maxSV/minSV

endif else begin

if (dim[0] ne dim[1]) then \$
MESSAGE, 'Input must be a square matrix.'
InverseA = INVERT(Array, Status, DOUBLE = dbl)

; Inversion failed?
if (Status eq 1) then return, -1

; Valid inverse.
norm1 = NORM(Array, DOUBLE=dbl, LNORM=lnorm)
norm2 = NORM(InverseA, DOUBLE=dbl, LNORM=lnorm)
return, norm1*norm2

endelse

end
```