/*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 "dcov3.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ dcov3(
double *v,
long mdim,
long n,
double sing[],
double var,
double work[],
long *ierr)
{
#define V(I_,J_)	(*(v+(I_)*(mdim)+(J_)))
	long int i, j;
	double stddev;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Sing = &sing[0] - 1;
	double *const Work = &work[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.
	 *>> 1996-03-30 DCOV3  Krogh   Added external statement.
	 *>> 1994-10-20 DCOV3  Krogh  Changes to use M77CON
	 *>> 1987-11-24 DCOV3  Lawson  Initial code.
	 *--D replaces "?": ?COV3, ?COPY, ?DOT, ?SCAL
	 *     Computes the covariance matrix for the solution vector of a
	 *     least-squares problem, Ax ~=~ b.  Assumes quantities are
	 *     available that have been computed by the singular value
	 *     decomposition subroutine, _SVDRS.
	 *     The covariance matrix is given by
	 *           VAR * (V*Pseudoinverse(S)) * Transp(V*Pseudoinverse(S))
	 *     ------------------------------------------------------------------
	 *                 Subroutine Arguments
	 *
	 *     V(,) [inout] Array containing the NxN orthogonal V matrix of the
	 *                  singular value decomposition of the matrix, A.
	 *                  On return contains the NxN symmetric
	 *                  covariance matrix.
	 *     MDIM   [in]  First dimension of the array, V(,).
	 *                  Require MDIM .ge. N.
	 *     N      [in]  Order of the matrix V contained in the array V(,).
	 *     SING() [in]  Singular values of the matrix, A.
	 *     VAR    [in]  Estimate of variance of error in the right-side
	 *                  vector, b, of the least-squares problem, Ax ~=~ b.
	 *     WORK() [scratch]  Work space of size N.
	 *     IERR   [out] Set to 0 if ok.  Set to J > 0 if the Jth singular
	 *                  value is zero.  In this latter case
	 *                  the covariance matrix cannot be computed and the
	 *                  contents of V(,) on return will be meaningless.
	 *     ------------------------------------------------------------------
	 *
	 *     May, 1987, C. L. Lawson & S. Y. Chiu, JPL.
	 *     Programmed in Fortran 77 for use in the JPL MATH77 library.
	 *     Prefixing subprogram names with S or D for s.p. or d.p. versions.
	 *     Using BLAS subprograms DCOPY, DDOT, & DSCAL, and
	 *     MATH77 error processing subr., IERM1, which uses IERMV1, ERMSG,
	 *     ERFIN, and AMACH.
	 *     ------------------------------------------------------------------ */
 
	/*     ------------------------------------------------------------------ */
	stddev = sqrt( var );
 
	/*              For J = 1, ...,N,  multiply col J of V by STDDEV/SING(J)
	 * */
	for (j = 1; j <= n; j++)
	{
		if (Sing[j] == ZERO)
		{
			ierm1( "DCOV3", 1, 0, "Jth singular value is zero", "J"
			 , j, '.' );
			*ierr = j;
			return;
		}
		dscal( n, stddev/Sing[j], &V(j - 1,0), 1 );
	}
 
	for (i = 1; i <= n; i++)
	{
		dcopy( n, &V(0,i - 1), mdim, work, 1 );
		for (j = i; j <= n; j++)
		{
			V(j - 1,i - 1) = ddot( n, work, 1, &V(0,j - 1), mdim );
		}
		for (j = 1; j <= (i - 1); j++)
		{
			V(j - 1,i - 1) = V(i - 1,j - 1);
		}
	}
	*ierr = 0;
	return;
#undef	V
} /* end of function */