/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:52 */
/*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 "dchol.h"
#include <float.h>
#include <stdlib.h>
void /*FUNCTION*/ dchol(
double *a,
long nda,
long n,
double b[],
double *s,
double tol,
long *ierr)
{
#define A(I_,J_)	(*(a+(I_)*(nda)+(J_)))
	long int _l0, i, j, k;
	double g, gmin, toli, tsq;
	static double tolm = 0.e0;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const B = &b[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.
	 *>> 2006-04-10 DCHOL Krogh  Fix for out of range unused subscript.
	 *>> 2001-10-25 DCHOL Snyder & Krogh Added c$OMP PARALLEL DO, etc.
	 *>> 2001-05-25 DCHOL Snyder & Krogh Compute inner products using DDOT
	 *>> 1996-03-30 DCHOL Krogh  Added external statement.
	 *>> 1994-10-20 DCHOL Krogh  Changes to use M77CON
	 *>> 1991-06-26 DCHOL Krogh  Initial MATH77 version. */
	/*     C. L. Lawson, JPL, 1970 April 30
	 *     Extensively modified by F. T. Krogh, JPL, 1991, September 23. */
	/* Solution of positive definite system using Cholesky decomposition. */
	/* Suppose the overdetermined system to be solved in the sense of least
	 * squares is   C*X = D */
	/* Suppose normal equations are formed and stored as follows. */
	/*                    A = (C**T)*C
	 *                    B = (C**T)*D
	 *                    S = (D**T)*D */
	/* Given B and the upper triangle of A this subroutine computes the
	 * solution vector X, storing it in place of B. */
	/* On input S should either = (D**T)D, or be zero.  If the former, then
	 * on return it is the euclidean norm of (CX - D).  Else the zero value
	 * is returned. */
	/* The upper triangle of A is replaced by the upper triangular Cholesky
	 * matrix, R, satisfying  (R**T)*R = A. */
	/* The value of TOL should be 10**(-k-1) where k is the minimum number of
	 * significant digits in A and B.   If TOL is < the relative precision in
	 * the computer's arithmetic, it is effectively replaced by that relative
	 * precision.  If some R(I,I)**2 is < TOL * (corresponding diagonal of
	 * A), then IERR is set to the value of I for the algebraically smallest
	 * value of R(I,I).  If R(I,I)**2 is .le. zero, then everything in the
	 * I-th row of R is set to zero, as is the I-th component of the
	 * solution, and IERR is replaced by -IERR to flag that this occured. */
	/*--D replaces "?": ?CHOL, ?DOT */
	/* ********************* Variable definitions *************************** */
	/* A      Input matrix, only upper triangle is used.  Replace by factor.
	 * B      Input right hand side, output solution.
	 * D1MACH Library routine to get characteristics of floating point
	 *  numbers.  D1MACH(4) = smallest number, x, such that 1.0 + x .ne. 1.0.
	 * DDOT   Compute dot products.
	 * G      Temporary for accumulating sums.
	 * GMIN   Smallest value seen for G when working on a diagonal entry.
	 * I      Index used to access matrix entries.
	 * IERR   Used to return status, see above.
	 * J      Index used to access matrix entries.
	 * K      Index used to access matrix entries.
	 * N      Order of matrix.
	 * NDA    Declared first dimension for A.
	 * S      Formal argument, see above.
	 * TOL    Input value for the tolerance, see above.
	 * TOLI   Internal value for the tolerance.
	 * TOLM   Machine tolerance, lower limit for TOLI.
	 * TSQ    TOLI**2. */
	/* ******************** Variable declarations *************************** */
	/* ******************** Start of executable code ************************ */
	if (n <= 0)
		return;
	if (tolm <= 0.e0)
		tolm = DBL_EPSILON;
	toli = fmax( tol, tolm );
	tsq = SQ(toli);
	gmin = A(0,0);
	*ierr = 0;
	for (i = 1; i <= n; i++)
	{
		g = A(i - 1,i - 1) - ddot( i - 1, &A(i - 1,0), 1, &A(i - 1,0),
		 1 );
		if (g < tsq*fabs( A(i - 1,i - 1) ))
		{
			if (g <= gmin)
			{
				gmin = g;
				*ierr = i;
				if (g <= 0.e0)
					*ierr = -i;
			}
			if (g <= 0.e0)
			{
				for (k = i; k <= n; k++)
				{
					A(k - 1,i - 1) = 0.e0;
				}
				B[i] = 0.e0;
				goto L_50;
			}
		}
		A(i - 1,i - 1) = sqrt( g );
 
		/*$OMP PARALLEL DO */
		for (j = i + 1; j <= n; j++)
		{
			A(j - 1,i - 1) = (A(j - 1,i - 1) - ddot( i - 1, &A(i - 1,0),
			 1, &A(j - 1,0), 1 ))/A(i - 1,i - 1);
		}
		/*$OMP END PARALLEL DO
		 *                        Solve next row of first lower triangular system */
		B[i] = (B[i] - ddot( i - 1, &A(i - 1,0), 1, b, 1 ))/A(i - 1,i - 1);
L_50:
		;
	}
	/*                        Get the solution norm */
	if (*s > 0.e0)
		*s = sqrt( fmax( 0.e0, *s - ddot( n, b, 1, b, 1 ) ) );
	/*                        Solve the second lower triangular system. */
	if (A(n - 1,n - 1) > 0.e0)
		B[n] /= A(n - 1,n - 1);
	for (i = n - 1; i >= 1; i--)
	{
		if (A(i - 1,i - 1) > 0.e0)
		{
			B[i] = (B[i] - ddot( n - i, &A(i,i - 1), nda, &B[i + 1],
			 1 ))/A(i - 1,i - 1);
		}
	}
	return;
#undef	A
} /* end of function */