/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:44 */
/*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 "durev.h"
#include <stdlib.h>
	/* COMMON translations */
struct t_ucom1 {
	long int n, m1, m2;
	}	ucom1;
	/* end of COMMON translations */
void /*FUNCTION*/ durev(
double *xt,
double *tx,
long ndim,
double *rcond,
long iwork[],
double *work)
{
#define XT(I_,J_)	(*(xt+(I_)*(ndim)+(J_)))
#define TX(I_,J_)	(*(tx+(I_)*(ndim)+(J_)))
#define WORK(I_,J_,K_)	(*(work+(I_)*(ucom1.n)*(ucom1.n)+(J_)*(ucom1.n)+\
	 (K_)))
	long int i, it, itxx, ix, ixtt, j, k;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	long *const Iwork = &iwork[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 DUREV Krogh  Added external statement.
	 *>> 1994-10-20 DUREV Krogh  Changes to use M77CON
	 *>> 1994-08-04 DUREV CLL Changed name of common from UCOM to UCOM1.
	 *>> 1992-02-17 CLL
	 *>> 1990-12-12 CLL DUREV  Initial code.
	 *     This subr does what is sometimes called series reversion.
	 *     Regarding N variables x sub j as functions of
	 *     N variables t sub i, and given a set of values of the t's and
	 *     the values of the x's and the 1st and 2nd partial derivatives of
	 *     the x's with respect to the t's evaluated at this set of t values,
	 *     this subr computes values at this point of the 1st and 2nd partial
	 *     derivatives of the t's with respect to the x's.
	 *     It is required that the matrix of 1st partials of the x's with
	 *     respect to the t's must be nonsingular.
	 *     ------------------------------------------------------------------
	 *                   Subroutine Arguments
	 *
	 *  XT(,) [in float]  Array containing values of N variables x sub 1
	 *     through x sub N, along with their partial derivatives of orders
	 *     1 and 2 with respect to N variables t sub 1 through t sub N, all
	 *     evaluated at the set of t values given in TX().
	 *     Data for x sub i is in XT(1:(N+2)*(N+1)/2),i).
	 *
	 *  TX(,) [inout float]  On entry TX(1,i), i = 1,..., N, must contain
	 *     values of the N variables t sub 1 thru t sub N.
	 *     Other values in the array TX() on entry are irrelevant.
	 *     This subr will compute 1st and 2nd partials of the t's with
	 *     respect to the x's and store these results in
	 *     TX(2:(N+2)*(N+1)/2),1:N).
	 *
	 *  NDIM  [in, integer]  Leading dimension for the arrays XT() and TX().
	 *     Require NDIM .ge. (N+2)*(N+1)/2.
	 *
	 *  RCOND  [out float]  Estimate of the reciprocal condition number of
	 *     the matrix of 1st partials of x's with respect to the t's.
	 *     RCOND will satisfy 0.0 .le. RCOND .le. 1.0.
	 *     Values near 1.0 indicate a well conditioned matrix, small values
	 *     indicate poor conditioning, zero indicates a singular matrix.  In
	 *     this latter case the subr will return without computing partials
	 *     of the t's with respect to the x's.
	 *
	 *  IWORK()  [scratch, integer]  Integer work space for this subroutine.
	 *
	 *  WORK(N,N,3) [scratch, float]  Floating-point work space for this
	 *     subroutine.  We use this as 3 NxN arrays.  In the calls to DGECO
	 *     and DGEI the last argument needs only N scratch locations.
	 *     ------------------------------------------------------------------
	 *              N, M1, M2 in common /UCOM1/
	 *
	 *  N [in]  Number of components in the (conceptual) x and t vectors.
	 *  M1, M2 [in]  We assume 0 .le. M1 .le. M2 .le. 2.  In other subrs of
	 *     the [D/S]UCOMP package, M1 and M2 select computation of partial
	 *     derivatives of orders M1 through M2, assuming all needed partial
	 *     derivatives of orders less than M1 are available.
	 *     This subr differs from others in the [D/S]UCOMP package in that
	 *     derivs of order 0, i.e., the values of t sub i for i = 1, ..., N,
	 *     associated with x sub j for j = 1, ..., N, must always be input to
	 *     this subr.  Thus this subr treats M1 = 0 like M1 = 1.
	 *     If M2 .eq. 0 this subr does nothing.
	 *     If M2 .eq. 1 this subr computes first partials of t w.r.t. x.
	 *     If M1 .le. 1 and M2 .eq. 2 this subr computes first and second
	 *     partials of t w.r.t. x.
	 *     If M1 .eq. 2 and M2 .eq. 2 this subr assumes first partials of t
	 *     w.r.t. x are available and computes second partials of t w.r.t. x.
	 *     ------------------------------------------------------------------
	 *--D replaces "?": ?UREV, ?DOT, ?GECO, ?GEI
	 *     Also uses ERMSG
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	if (ucom1.m2 == 0)
		return;
	if (ucom1.m1 <= 1)
	{
		/*          Copy the matrix of 1st partials of x w.r.t. t to WORK2(*,*,1)
		 * */
		for (ix = 1; ix <= ucom1.n; ix++)
		{
			for (it = 1; it <= ucom1.n; it++)
			{
				WORK(0,it - 1,ix - 1) = XT(ix - 1,it);
			}
		}
 
		/*        Compute an LU factorization and RCOND for this Jacobian matrix.
		 * */
		dgeco( &WORK(0,0,0), ucom1.n, ucom1.n, iwork, rcond, &WORK(2,0,0) );
		if (*rcond == 0.0e0)
		{
			ermsg( "DUREV", 6, 0, "Singular Jacobian matrix.", '.' );
			return;
		}
 
		/*        Compute the inverse Jacobian matrix in WORK(*,*,1) and copy
		 *        its elements into TX().  These elements are the 1st partials of
		 *        the t's w.r.t. the x's.
		 * */
		dgei( &WORK(0,0,0), ucom1.n, ucom1.n, iwork, &WORK(2,0,0) );
 
		for (it = 1; it <= ucom1.n; it++)
		{
			for (ix = 1; ix <= ucom1.n; ix++)
			{
				TX(it - 1,ix) = WORK(0,ix - 1,it - 1);
			}
		}
	}
	else
	{
 
		/*          Copy 1st partials of t's w.r.t. x's from TX() to WORK(*,*,1).
		 * */
		for (it = 1; it <= ucom1.n; it++)
		{
			for (ix = 1; ix <= ucom1.n; ix++)
			{
				WORK(0,ix - 1,it - 1) = TX(it - 1,ix);
			}
		}
	}
	/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
	if (ucom1.m2 < 2)
		return;
	/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	 *                 Main loop on K to compute 2nd partials
	 *                 of t sub K w.r.t. the x's.
	 * */
	for (k = 1; k <= ucom1.n; k++)
	{
 
		/*           Compute weighted sum of the Hessian matrices of the x's
		 *           using the 1st partial of t sub K w.r.t. x sub i as the
		 *           weight for the Hessian matrix of x sub i.
		 *           Store the resulting symmetric matrix in WORK(*,*,2).
		 * */
		ixtt = 1 + ucom1.n;
		for (i = 1; i <= ucom1.n; i++)
		{
			for (j = 1; j <= i; j++)
			{
				ixtt += 1;
				WORK(1,j - 1,i - 1) = ddot( ucom1.n, &WORK(0,0,k - 1),
				 ucom1.n, &XT(0,ixtt - 1), ndim );
				WORK(1,i - 1,j - 1) = WORK(1,j - 1,i - 1);
			}
		}
 
		/*           Multiply the matrix in WORK(*,*,2) by WORK(*,*,1) from the
		 *           right and by the transpose of WORK(*,*,1) from the left and
		 *           invert the sign.  Put the result of the first multiplication
		 *           into WORK(*,*,3) and the lower triangle of the (symmetric)
		 *           final result into TX(*,K).
		 * */
		for (i = 1; i <= ucom1.n; i++)
		{
			for (j = 1; j <= ucom1.n; j++)
			{
				WORK(2,j - 1,i - 1) = ddot( ucom1.n, &WORK(1,0,i - 1),
				 ucom1.n, &WORK(0,j - 1,0), 1 );
			}
		}
 
		itxx = 1 + ucom1.n;
		for (i = 1; i <= ucom1.n; i++)
		{
			for (j = 1; j <= i; j++)
			{
				itxx += 1;
				TX(k - 1,itxx - 1) = -ddot( ucom1.n, &WORK(0,i - 1,0),
				 1, &WORK(2,j - 1,0), 1 );
			}
		}
	}
	return;
#undef	WORK
#undef	XT
#undef	TX
} /* end of function */