;$Id: cond.pro,v 1.11 1997/01/15 03:11:50 ali Exp $
;
; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
;       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.
;
; 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) by explicitly computing, 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 nr_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.PRO results in floating-
;       point errors.
;
; REFERENCE:
;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
;       Erwin Kreyszig
;       ISBN 0-471-55380-8
;
; 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. 
;-

FUNCTION Cond, A, Double = Double

  ON_ERROR, 2  ;Return to caller if error occurs.

  Dimension = SIZE(A)
  if Dimension[0] ne 2 then MESSAGE, $
    "Input array must be two-dimensional."

  if Dimension[1] ne Dimension[2] then MESSAGE, $
    "Input array must be square."

  if N_ELEMENTS(Double) eq 0 then $
    Double = (Dimension[Dimension[0]+1] eq 5) or $
             (Dimension[Dimension[0]+1] eq 9)

  InverseA = INVERT(A, Status, Double = Double)
  if Status eq 0 then RETURN, $ ;Valid inverse.
    NORM(A, Double = Double) * NORM(InverseA, Double = Double) $
  else RETURN, -1 ;Inversion failed.

END