/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:18 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_snqsol s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include "p_snqsol.h"
/*     program DRSNQSOL
 *>> 1996-06-21 DRSNQSOL Krogh  Changes for C conversion.
 *>> 1994-11-02 DRSNQSOL Krogh  Changes to use M77CON
 *>> 1992-04-15 DRSNQSOL CLL.
 *>> 1992-01-14 CLL.
 *     Demo driver for SNQSOL.  Also using SCKDER.
 *     Expected solution vector:  0.9000518     1.0001835     1.0945009
 *     ------------------------------------------------------------------
 *--S replaces "?": DR?NQSOL, ?NQSOL, ?NRM2, ?CKDER, ?NQFJ
 *     ------------------------------------------------------------------ */
		/* PARAMETER translations */
#define	LWA	(3 + (15*NMAX + 3*NMAX*NMAX)/2)
#define	NMAX	3
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int _l0, i, imax, iopt[5], j, jmax, m, mode;
	float fjac[NMAX][NMAX], fnorm, fvec[NMAX], test[NMAX][NMAX],
	 tol, tstmax, wa[LWA], x[NMAX];
	static long n = NMAX;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Fvec = &fvec[0] - 1;
	long *const Iopt = &iopt[0] - 1;
	float *const Wa = &wa[0] - 1;
	float *const X = &x[0] - 1;
		/* end of OFFSET VECTORS */
 
	/*     ------------------------------------------------------------------ */
	Iopt[4] = 0;
	tol = sqrtf( FLT_EPSILON );
	X[1] = 3.0e0;
	X[2] = 3.0e0;
	X[3] = 3.0e0;
	printf("Program DRSNQSOL.  Demo driver for SNQSOL.  Also using SCKDER.\n");
 
	/*     ------------------------------------------------------------------
	 *               Using SCKDER to check derivative computation.
	 *     ------------------------------------------------------------------ */
	printf("\n Using SCKDER to check derivative computation.\n");
	m = n;
	snqfj( n, x, fvec, (float*)fjac, ADR(_l0,2) );
	mode = 1;
L_10:
	;
	sckder( &mode, m, n, x, fvec, (float*)fjac, NMAX, (float*)test,
	 &imax, &jmax, &tstmax );
	if (mode == 2)
	{
		snqfj( n, x, fvec, (float*)fjac, ADR(_l0,1) );
		goto L_10;
	}
	snqfj( n, x, fvec, (float*)fjac, ADR(_l0,1) );
	printf("\n           X(J) =");
	for (j = 1; j <= n; j++)
	{
		printf("%11.3g", X[j]);
	}
	printf("\n");
	printf("\n   I    FVEC(I)   .....................FJAC(I,J)........................\n\n");
	for (i = 1; i <= m; i++)
	{
		printf(" %3ld %11.3g", i, Fvec[i]);
		printf(" ");
		for (j = 1; j <= n; j++)
		{
			printf("%11.3g", fjac[j - 1][i - 1]);
		}
		printf("\n");
	}
	printf("\n TEST(,):\n\n");
	for (i = 1; i <= m; i++)
	{
		printf(" %3ld", i);
		printf("             ");
		for (j = 1; j <= n; j++)
		{
			printf("%11.3g", test[j - 1][i - 1]);
		}
		printf("\n");
	}
	printf("\n IMAX =%3ld,    JMAX =%3ld,    TSTMAX =%11.3g\n", imax, jmax, tstmax);
 
	/*     ------------------------------------------------------------------
	 *           Using SNQSOL to solve system of nonlinear equations.
	 *     ------------------------------------------------------------------ */
	printf("\n Using SNQSOL to solve system of nonlinear equations.\n");
 
	snqsol( snqfj, n, x, fvec, tol, iopt, wa, LWA );
 
	fnorm = snrm2( n, fvec, 1 );
	/*++ CODE for ~.C. is inactive
	 *      print'('' Termination status:  '',i6/'' NFEV, NJEV:          '',
	 *     * 2i6/ '' Final residual norm: '',g14.3/'' Final X():           ''
	 *     * /(8x,4f14.7))',  IOPT(1), IOPT(2), IOPT(3),
	 *     * FNORM, (X(J), J = 1, N)
	 *++ CODE for .C. is active */
   printf(" Termination status:  %6ld\n NFEV, NJEV:          "
      "%6ld%6ld\n Final residual norm: %14.3g\n Final X():           ", Iopt[1], Iopt[2],
      Iopt[3], fnorm);
   for (j = 0; j < n; j+=4){
      printf("\n        ");
      for (i = j; i < (j < n - 3 ? j + 4 : n); i++)
         printf("%14.7f", x[i]);}
   printf("\n");
	exit(0);
	/*++ END */
} /* end of function */
/*     ================================================================== */
void /*FUNCTION*/ snqfj(
long n,
float x[],
float fvec[],
float *fjac,
long *iflag)
{
#define FJAC(I_,J_)	(*(fjac+(I_)*(n)+(J_)))
	long int i;
	static float c1[3]={-1.0e0,2.0e0,2.0e0};
	static float c2[3]={2.0e0,-1.0e0,2.0e0};
	static float c3[3]={2.0e0,2.0e0,-1.0e0};
	static float term[3]={5.01e0,5.85e0,8.88e0};
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const C1 = &c1[0] - 1;
	float *const C2 = &c2[0] - 1;
	float *const C3 = &c3[0] - 1;
	float *const Fvec = &fvec[0] - 1;
	float *const Term = &term[0] - 1;
	float *const X = &x[0] - 1;
		/* end of OFFSET VECTORS */
 
	/*>> 1992-01-14 CLL.
	 *     Sample 3-dimensional function of 3 variables for demo of solution
	 *     of a system of nonlinear equations.
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	if (*iflag == 1)
	{
		/*                                        Compute function vector. */
		for (i = 1; i <= n; i++)
		{
			Fvec[i] = expf( C1[i]*X[1] ) + sinhf( C2[i]*X[2] ) + tanhf( C3[i]*
			 X[3] ) - Term[i];
		}
	}
	else if (*iflag == 2)
	{
		/*                                        Compute Jacobian matrix. */
		for (i = 1; i <= n; i++)
		{
			FJAC(0,i - 1) = expf( C1[i]*X[1] )*C1[i];
			FJAC(1,i - 1) = coshf( C2[i]*X[2] )*C2[i];
			FJAC(2,i - 1) = powif(1.0e0/coshf( C3[i]*X[3] ),2)*C3[i];
		}
	}
	return;
#undef	FJAC
} /* end of function */