;+
; NAME:
;   MCHOLDC
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Modified Cholesky Factorization of a Symmetric Matrix
;
; MAJOR TOPICS:
;   Linear Systems
;
; CALLING SEQUENCE:
;   MCHOLDC, A, D, E [, /OUTFULL, /SPARSE, /PIVOT, TAU=TAU, $
;                      PERMUTE=PERMUTE, INVPERMUTE=INVPERMUTE ]
;
; DESCRIPTION:
;
;  Given a symmetric matrix A, the MCHOLDC procedure computes the
;  factorization:
;
;     A + E   =   TRANSPOSE(U) ## D ## U
;
;  where A is the original matrix (optionally permuted if the PIVOT
;  keyword is set), U is an upper triangular matrix, D is a diagonal
;  matrix, and E is a diagonal error matrix.
;
;  The standard Cholesky factorization is only defined for a positive
;  definite symmetric matrix.  If the input matrix is positive
;  definite then the error term E will be zero upon output.  The user
;  may in fact test the positive-definiteness of their matrix by
;  factoring it and testing that all terms in E are zero.
;
;  If A is *not* positive definite, then the standard Cholesky
;  factorization is undefined.  In that case we adopt the "modified"
;  factorization strategy of Gill, Murray and Wright (p. 108), which
;  involves adding a diagonal error term to A in order to enforce
;  positive-definiteness.  The approach is optimal in the sense that
;  it attempts to minimize E, and thus disturbing A as little as
;  possible.  For optimization problems, this approximate
;  factorization can be used to find a direction of descent even when
;  the curvature is not positive definite.
;
;  The upper triangle of A is modified in place.  By default, the
;  lower triangle is left unchanged, and the matrices D and E are
;  actually returned as vectors containing only the diagonal terms.
;  However, if the keyword OUTFULL is set then full matrices are
;  returned.  This is useful when matrix multiplication will be
;  performed at the next step.
;
;  The modified Cholesky factorization is most stable when pivoting is
;  performed.  If the keyword PIVOT is set, then pivoting is performed
;  to place the diagonal terms with the largest amplitude in the next
;  row.  The permutation vectors returned in PERMUTE and INVPERMUTE
;  can be used to apply and reverse the pivoting.
;    [ i.e.,  (U(PP,*))(*,PP) applies the permutation and
;             (U(IPP,*))(*,IPP) reverses it, where PP and IPP are the
;             permutation and inverse permutation vectors. ]
;
;  If the matrix to be factored is very sparse, then setting the
;  SPARSE keyword may improve the speed of the computations.  SPARSE
;  is more costly on a dense matrix, but only grows as N^2, where as
;  the standard computation grows as N^3, where N is the rank of the
;  matrix.
;
;  If the CHOLSOL keyword is set, then the output is slightly
;  modified.  The returned matrix A that is returned is structured so
;  that it is compatible with the CHOLSOL built-in IDL routine.  This
;  involves converting A to being upper to lower triangular, and
;  multiplying by SQRT(D).  Users must be sure to check that all
;  elements of E are zero before using CHOLSOL.
;
; PARAMETERS:
;
;   A - upon input, a symmetric NxN matrix to be factored.
;       Upon output, the upper triangle of the matrix is modified to
;       contain the factorization.
;
;   D - upon output, the diagonal matrix D.
;
;   E - upon output, the diagonal error matrix E.
;
; KEYWORD PARAMETERS:
;
;   OUTFULL - if set, then A, D and E will be modified to be full IDL
;             matrices than can be matrix-multiplied.  By default,
;             only the upper triangle of A is modified, and D and E
;             are returned as vectors.
;
;   PIVOT - if set, then diagonal elements of A are pivoted into place
;           and operated on, in decrease order of their amplitude.
;           The permutation vectors are returned in the PERMUTE and
;           INVPERMUTE keywords.
;
;   PERMUTE - upon return, the permutation vector which converts a
;             vector into permuted form.
;
;   INVPERMUTE - upon return, the inverse permutation vector which
;                converts a vector from permuted form back into
;                standard form.
;
;   SPARSE - if set, then operations optimized for sparse matrices are
;            employed.  For large but very sparse matrices, this can
;            save a significant amount of computation time.
;
;   CHOLSOL - if set, then A and D are returned, suitable for input to
;             the built-in IDL routine CHOLSOL.  CHOLSOL is mutually
;             exclusive with the FULL keyword.
;
;   TAU - if set, then use the Tau factor as described in the
;         "unconventional" modified Cholesky factorization, as
;         described by Xie & Schlick.
;         Default: the unconventional technique is not used.
;
; EXAMPLE:
;
;   Example 1
;   ---------
;   a = randomn(seed, 5,5)    ;; Generate a random matrix
;   a = 0.5*(transpose(a)+a)  ;; Symmetrize it
;
;   a1 = a                    ;; Make a copy
;   mcholdc, a1, d, e, /full  ;; Factorize it
;   print, max(abs(e))        ;; This matrix is not positive definite
;
;   diff = transpose(a1) ## d ## a1 - e - a
;                             ;; Test the factorization by inverting
;                             ;; it and subtracting A
;   print, max(abs(diff))     ;; Differences are small
;
;   Example 2
;   ---------
;   Solving a problem with MCHOLDC and CHOLSOL
;
;   a = [[6E,15,55],[15E,55,225],[55E,225,979]]
;   b = [9.5E,50,237]
;
;   mcholdc, a, d, e, /cholsol  ;; Factorize matrix, compatible w/ CHOLSOL
;   if total(abs(e)) NE 0 then $
;      message, 'ERROR: Matrix A is not positive definite'
;
;   x = cholsol(a, d, b)        ;; Solve with CHOLSOL
;   print, x
;        -0.500001    -0.999999     0.500000
;   which is within 1e-6 of the true solution.
;
;
; REFERENCES:
;
;   Gill, P. E., Murray, W., & Wright, M. H. 1981
;     *Practical Optimization*, Academic Press
;   Schlick, T. & Fogelson, A., "TNPACK - A Truncated Newton
;     Minimization Package for Large- Scale Problems: I. Algorithm and
;     Usage," 1992, ACM TOMS, v. 18, p. 46-70.  (Alg. 702)
;   Xie, D. & Schlick, T., "Remark on Algorithm 702 - The Updated
;     Truncated Newton Minimization Package," 1999, ACM TOMS, v. 25,
;     p. 108-122.
;
; MODIFICATION HISTORY:
;   Written, CM, Apr 2001
;   Added CHOLSOL keyword, CM, 15 Feb 2002
;   Fix bug when computing final correction factor (thanks to Aaron
;     Adcock), CM, 13 Nov 2010
;
;  $Id: mcholdc.pro,v 1.6 2010/11/13 09:45:55 cmarkwar Exp $
;
;-
; Copyright (C) 2001, 2002, 2010, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-