/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:40 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "dcov2.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ONE	1.0e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ dcov2(
double *a,
long mdim,
long n,
long ip[],
double var,
long *ierr)
{
#define A(I_,J_)	(*(a+(I_)*(mdim)+(J_)))
	long int i, ip1, j, k, kp1;
	double tmp;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	long *const Ip = &ip[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 2001-02-13 DCOV2  Krogh   Fixed a comment.
	 *>> 1996-03-30 DCOV2  Krogh   Added external statement.
	 *>> 1994-11-11 DCOV2  Krogh   Declared all vars.
	 *>> 1994-10-20 DCOV2  Krogh  Changes to use M77CON
	 *>> 1989-10-20 DCOV2  CLL
	 *>> 1987-11-24 DCOV2  Lawson  Initial code.
	 *--D replaces "?": ?COV2, ?DOT, ?SWAP
	 *     Computes upper triangle of covariance matrix,
	 *     beginning with the triangular matrix, and permutation record, and
	 *     data error variance computed by _HFTI.
	 *     Thus, given a matrix, A, represented by the upper triangular
	 *     matrix in A(), and the permutation record in IP(), and the data
	 *     error variance, VAR, we compute the upper triangle of the
	 *     symmetric matrix, C = VAR * ((A**t)*A)**(-1).
	 *     Adapted from PROG2 in L & H book.
	 *     ------------------------------------------------------------------
	 *     Subprograms referenced directly: DDOT, DSWAP, IERM1
	 *     Other subprograms needed: ERMSG, IERV1, ERFIN, AMACH
	 *     ------------------------------------------------------------------
	 *                 Subroutine Arguments
	 *
	 *     A(,) [inout] On entry, contains the upper triangular matrix, A,
	 *                  in standard, not packed, storage.  This matrix could
	 *                  have been produced by _HFTI.  On return, contains the
	 *                  upper triangle of the symmetric unscaled covariance
	 *                  matrix.  Elements below the diagonal in A(,) will
	 *                  not be referenced.
	 *     MDIM [in]    First dimension of A(,).  Require MDIM .ge. N.
	 *     N [in]       Order of the matrix in A(,)
	 *     IP() [in]    Permutation record produced by _HFTI.
	 *     VAR  [in]    Estimate of variance of data error.
	 *     IERR [out]   Set to 0 if ok.  Set to J > 0 if the (J,J) element
	 *                  of the given matrix is zero.  In this latter case
	 *                  the covariance matrix cannot be computed and the
	 *                  contents of A(,) on return will be meaningless.
	 *     ------------------------------------------------------------------
	 *          This code was originally developed by Charles L. Lawson and
	 *     Richard J. Hanson at Jet Propulsion Laboratory in 1973.  The
	 *     original code was described and listed in the book,
	 *
	 *                  Solving Least Squares Problems
	 *                  C. L. Lawson and R. J. Hanson
	 *                  Prentice-Hall, 1974
	 *
	 *     Feb, 1985, C. L. Lawson & S. Y. Chan, JPL.  Adapted code from the
	 *     Lawson & Hanson book to Fortran 77 for use in the JPL MATH77
	 *     library.
	 *     Prefixing subprogram names with S or D for s.p. or d.p. versions.
	 *     Using generic names for intrinsic functions.
	 *     Adding calls to BLAS and MATH77 error processing subrs in some
	 *     program units.
	 *     1989-10-20 CLL Moved integer declaration earlier to avoid warning
	 *     msg from Cray compiler.
	 *     ------------------------------------------------------------------ */
 
	/*     ------------------------------------------------------------------
	 *     Replace upper triangular matrix U by its inverse, Inv(U)
	 * */
	for (j = 1; j <= n; j++)
	{
		if (A(j - 1,j - 1) == ZERO)
		{
			ierm1( "DCOV2", 1, 0, "Jth diagonal elt is zero", "J",
			 j, '.' );
			*ierr = j;
			return;
		}
		A(j - 1,j - 1) = ONE/A(j - 1,j - 1);
	}
 
	for (i = 1; i <= (n - 1); i++)
	{
		for (j = i + 1; j <= n; j++)
		{
			A(j - 1,i - 1) = -A(j - 1,j - 1)*ddot( j - i, &A(i - 1,i - 1),
			 mdim, &A(j - 1,i - 1), 1 );
		}
	}
 
	/*     Replace Inv(U) by upper triangle of Inv(u) * Trans(Inv(U))
	 *     multiplied by VAR.
	 * */
	for (i = 1; i <= n; i++)
	{
		for (j = i; j <= n; j++)
		{
			A(j - 1,i - 1) = var*ddot( n - j + 1, &A(j - 1,i - 1),
			 mdim, &A(j - 1,j - 1), mdim );
		}
	}
	/*                                 Permute rows & columns */
	for (i = n - 1; i >= 1; i--)
	{
		if (Ip[i] != i)
		{
			k = Ip[i];
			tmp = A(i - 1,i - 1);
			A(i - 1,i - 1) = A(k - 1,k - 1);
			A(k - 1,k - 1) = tmp;
			if (i > 1)
			{
				dswap( i - 1, &A(i - 1,0), 1, &A(k - 1,0), 1 );
			}
			ip1 = i + 1;
			if (ip1 < k)
			{
				dswap( k - i - 1, &A(ip1 - 1,i - 1), mdim, &A(k - 1,ip1 - 1),
				 1 );
			}
			kp1 = k + 1;
			if (kp1 <= n)
			{
				dswap( n - k, &A(kp1 - 1,i - 1), mdim, &A(kp1 - 1,k - 1),
				 mdim );
			}
		}
	}
	*ierr = 0;
	return;
#undef	A
} /* end of function */