/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:04 */
/*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 "srangv.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ srangv(
float *a,
long ndim,
long n,
float u[],
float x[],
LOGICAL32 *havec,
long *ierr)
{
#define A(I_,J_)	(*(a+(I_)*(ndim)+(J_)))
	long int j, k, l;
	float sum;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const U = &u[0] - 1;
	float *const X = &x[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 SRANGV  Krogh   Added external statement.
	 *>> 1994-11-11 SRANGV Krogh   Declared all vars.
	 *>> 1994-10-20 SRANGV Krogh  Changes to use M77CON
	 *>> 1987-04-22 SRANGV Lawson  Initial code.
	 *--S replaces "?": ?RANGV, ?RANG
	 *
	 *     SRANGV generates an N-component random vector, X, from a multi-
	 *     variate normal distribution having a specified mean vector, U,
	 *     and covariance matrix, A.
	 *     ------------------------------------------------------------------
	 *                        Subroutine arguments
	 *
	 *     A(,)  [inout]  Array with first dimension NDIM, and second
	 *           dimension at least N.  When HAVEC is false, A(,) contains
	 *           the given N x N symmetric covariance matrix.  This
	 *           subroutine will replace A by its lower-triangular Cholesky
	 *           factor, C, and set HAVEC = .true.
	 *           When HAVEC is true A(,) is assumed to contain the Cholesky
	 *           factor, C.
	 *           Note:  Only the lower triangle of the array A(,) is used by
	 *           this subroutine.
	 *
	 *     NDIM  [in]  First dimension of the array A(,).
	 *           Require NDIM .ge. N.
	 *
	 *     N  [in]  Order of the covariance matrix A and dimension of the
	 *           vectors U and X.
	 *
	 *     U()  [in]   Contains the N-dimensional mean vector.
	 *
	 *     X()  [out]  On return will contain the N-dimensional generated
	 *                 random vector.
	 *
	 *     HAVEC    [inout]  See description above for A(,).
	 *
	 *     IERR  [out]    IERR is only referenced when this subr is entered
	 *           with HAVEC = false.
	 *           In that case IERR will be set to zero if the
	 *           Cholesky factorization is successful.  If the factorization
	 *           is unsuccessful IERR will be set to the index of the row at
	 *           which failure is detected.  In this latter case the results
	 *           returned in A(,) and X() will not be valid.
	 *     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	 *     Subprogram referenced: SRANG
	 *     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	 *     Programmed for JPL by Carl Pitts, Heliodyne Corp., April, 1969.
	 *     Adapted to Fortran 77 for the JPL MATH 77 library by C. L. Lawson
	 *     and S. Y. Chiu, JPL, April 1987.
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	if (!*havec)
	{
		/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		 *     BEGIN PROCEDURE TO CALCULATE C-MATRIX, NOTE THE C-MATRIX
		 *     overwrites the lower triangle of the A matrix.
		 * */
		for (j = 1; j <= n; j++)
		{
			for (k = 1; k <= j; k++)
			{
				sum = ZERO;
				for (l = 1; l <= (k - 1); l++)
				{
					sum += A(l - 1,j - 1)*A(l - 1,k - 1);
				}
				if (j != k)
				{
					A(k - 1,j - 1) = (A(k - 1,j - 1) - sum)/A(k - 1,k - 1);
				}
				else
				{
					sum = A(k - 1,k - 1) - sum;
					if (sum <= ZERO)
					{
						*ierr = k;
						/*                                                Error return. */
						return;
					}
					A(k - 1,k - 1) = sqrtf( sum );
				}
			}
		}
		*havec = TRUE;
		*ierr = 0;
	}
	/*                            END PROCEDURE TO CALCULATE C-MATRIX
	 *     ------------------------------------------------------------------
	 *                            BEGIN PROCEDURE TO CALCULATE RANDOM VECTOR
	 *
	 *     LOOP TO FILL Y-VECTOR WITH NORMALLY DISTRIBUTED RANDOM NUMBERS
	 *     WITH ZERO MEAN AND UNIT VARIANCE  NOTE THE Y-VECTOR IS STORED
	 *     TEMPORARILY IN THE X-VECTOR
	 * */
	for (j = 1; j <= n; j++)
	{
		X[j] = srang();
	}
 
	/*     LOOP TO CALCULATE X-VECTOR
	 * */
	for (j = n; j >= 1; j--)
	{
		sum = U[j];
		for (k = 1; k <= j; k++)
		{
			sum += A(k - 1,j - 1)*X[k];
		}
		X[j] = sum;
	}
	return;
#undef	A
} /* end of function */