/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */
/*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 "sckder.h"
#include <float.h>
#include <stdlib.h>
void /*FUNCTION*/ sckder(
long *mode,
long m,
long n,
float x[],
float fvec[],
float *fjac,
long ldfjac,
float *test,
long *imax,
long *jmax,
float *tstmax)
{
#define FJAC(I_,J_)	(*(fjac+(I_)*(ldfjac)+(J_)))
#define TEST(I_,J_)	(*(test+(I_)*(ldfjac)+(J_)))
	static LOGICAL32 plus;
	long int _l0, i;
	static long int j;
	float eps, fac, tstij;
	static float alpha, delx, small, xj;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Fvec = &fvec[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 SCKDER Krogh  Added external statement.
	 *>> 1994-10-20 SCKDER Krogh  Changes to use M77CON
	 *>> 1992-01-16 SCKDER CLL@JPL
	 *>> 1991-12-16 CLL@JPL
	 *  SCKDER checks a Jacobian matrix of first partial derivatives for
	 *  consistency with central finite differences of function values.
	 *
	 *     Usage:
	 *
	 *        Set X() to a trial N-vector.
	 *        Compute FJAC(,) as the MxN Jacobian matrix of partials of
	 *           FVEC with respect to X, evaluated at X().
	 *        MODE = 1
	 *     10 call SCKDER(MODE,...)
	 *        if(MODE .eq. 2) then
	 *           Compute FVEC() as an M-vector of function values
	 *              evaluated at X().
	 *           go to 10
	 *        endif
	 *        Here the process is completed.
	 *        Results are in IMAX, JMAX, TSTMAX, and TEST(,).
	 *     ------------------------------------------------------------------
	 *                         Subroutine Arguments
	 *
	 *  MODE [inout]  On initial entry set MODE = 1.
	 *       SCKDER will return a number of times with MODE = 2.  The calling
	 *       code should compute FVEC() as a function of X() and call SCKDER
	 *       again, not altering MODE.  When SCKDER is finished it returns
	 *       with MODE = 3.
	 *  M [in]  Number of terms in FVEC() and number of rows of data
	 *          in FJAC(,).
	 *  N [in]  Number of terms in X() and number of cloumns of data
	 *          in FJAC(,).
	 *  X() [inout]  Initially must contain the x-vector around which the
	 *          testing will be done.  Contains perturbed x-vectors on
	 *          returns with MODE = 2.  On final return with MODE = 3,
	 *          contains the original x-vector exactly restored.
	 *  FVEC() [in]  The user stores function values in FVEC().
	 *  FJAC(,) [in]  The user stores the Jacobian matrix in FVEC(,).
	 *  LDFJAC [in]  Declared first dimension of the arrays FJAC(,) and
	 *               TEST(,).  Require LDFJAC .ge. M.
	 *  TEST(,) [out]  Array with same dimensions as FJAC().
	 *     For each i = 1, ..., M, and j = 1, ..., N, SCKDER sets TEST(i,j)
	 *     to the signed difference: FJAC(i,j) minus a central finite-
	 *     difference approximation to the first partial derivative of
	 *     f sub i with respect to x sub j.
	 *  IMAX, JMAX [out]  I and J indices of the value in TEST(,) of
	 *               largest magnitude.
	 *  TSTMAX [out]  Magnitude of the value in TEST(,) of largest magnitude,
	 *               i.e., TSTMAX = abs(TEST(IMAX, JMAX))
	 *     ------------------------------------------------------------------
	 *  R1MACH(1) is the underflow limit.
	 *  R1MACH(4) is the machine precision.
	 *     ------------------------------------------------------------------
	 *--S replaces "?": ?CKDER
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	switch (*mode)
	{
		case 1: goto L_10;
		case 2: goto L_20;
	}
	/*                               ANSI Standard Fortran 77 comes here
	 *                               if MODE is not 1 or 2.
	 * */
	ierm1( "SCKDER", 1, 0, "Require MODE = 1 or 2.", "MODE", *mode,
	 '.' );
	return;
 
L_10:
	;
	eps = FLT_EPSILON;
	alpha = powf(3.0e0*eps,0.333333e0);
	small = fmaxf( 1.0e5*FLT_MIN/alpha, eps*eps );
	*tstmax = 0.0e0;
	*imax = 0;
	*jmax = 0;
	*mode = 2;
	j = 0;
L_15:
	;
	j += 1;
	xj = X[j];
	if (fabsf( xj ) > small)
	{
		delx = alpha*xj;
	}
	else if (fabsf( xj ) > 0.0e0)
	{
		delx = alpha*small;
	}
	else
	{
		delx = alpha;
	}
	X[j] = xj + delx;
	plus = TRUE;
	return;
 
	/*           Here we return to the user's code for computation of
	 *           FVEC() using X().  Execution will resume here.
	 * */
L_20:
	;
	if (plus)
	{
		/*                         Using col N of TEST() to save FVEC()
		 *                         evaluated using XJ + DELX. */
		for (i = 1; i <= m; i++)
		{
			TEST(n - 1,i - 1) = Fvec[i];
		}
		X[j] = xj - delx;
		plus = FALSE;
		return;
		/*                  Returning again to the user's code for another
		 *                  evaluation of FVEC() using the perturbed X(). */
	}
 
	fac = 0.5e0/delx;
	for (i = 1; i <= m; i++)
	{
		tstij = FJAC(j - 1,i - 1) - fac*(TEST(n - 1,i - 1) - Fvec[i]);
		TEST(j - 1,i - 1) = tstij;
		if (fabsf( tstij ) > *tstmax)
		{
			*tstmax = fabsf( tstij );
			*imax = i;
			*jmax = j;
		}
	}
	/*                                     Restore X(J) */
	X[j] = xj;
	if (j < n)
		goto L_15;
	*mode = 3;
	return;
#undef	TEST
#undef	FJAC
} /* end of function */