/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:20 */
/*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 "silupmd.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
		/* PARAMETER translations */
#define	LTXTAA	1
#define	LTXTAB	53
#define	LTXTAC	92
#define	LTXTAD	136
#define	LTXTAE	162
#define	LTXTAF	232
#define	LTXTAG	286
#define	LTXTAH	325
#define	LTXTAI	375
#define	LTXTAJ	408
#define	LTXTAK	455
#define	LTXTAL	508
#define	MAXDIM	10
#define	MEFVEC	61
#define	MEIVEC	57
#define	MEMDA1	27
#define	MENTXT	23
#define	MERET	51
#define	METDIG	22
#define	METEXT	53
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ silupmd(
long ndim,
float x[],
float y,
long ntab[],
float xt[],
float yt[],
long ndeg[],
long lup[],
long iopt[],
float eopt[])
{
	long int i, iiflg, in[MAXDIM], intchk[4], iyt, j, k, knt[MAXDIM],
	 l, lr[MAXDIM - 1], lt, lx[MAXDIM], ndimi, ntabm, ntabxt, nup[MAXDIM];
	float fdat[MAXDIM];
	static char mtxtaa[3][206]={"Inputs to SILUPM, with NDIM = $I, EOPT dim. of $I.$EIOPT($I) = $I, is not a valid option.$EIOPT($I) = $I, requests an error estimate.$EIOPT($I) = $I, specifies$EIOPT($I) = $I, Store derivatives at $I, Max.",
	 " total derivative of $I,$EIOPT($I) = $I, requests Abs. & Rel. errors of: $F $F$EIOPT($I) = $I, request accuracy of $F$EIOPT($I) = $I, skips data points with YT(?) = $F$ELUP($I) = $I; it should be < 4.$ERag",
	 "ged table badly specified in dimension $I.$EStart of ragged table, NTAB($I) = $I; it must be 0.$EMiddle of ragged table, NTAB($I) = $I; it must be >0.$EEnd of ragged table, NTAB($I) = $I; it must be -1.$E "};
	static char mtxtab[1][35]={"Extrpolation polynomial degrees:$B"};
	static char mtxtac[1][28]={"max derivatives per X(I):$B"};
	static char mtxtad[1][10]={"NDEG():$B"};
	static char mtxtae[1][9]={"LUP():$B"};
	static char mtxtaf[1][10]={"NTAB():$B"};
	static char mtxtag[1][8]={"X(): $B"};
	static char mtxtah[1][17]={"Dim. 1 XT's = $B"};
	static char mtxtai[1][23]={"Dim. 1 to $I XT's = $B"};
	static char mtxtaj[1][29]={"Dim. $I XT's = $F + k * %F$E"};
	static char mtxtak[1][18]={"Dim. $I XT's = $B"};
	static char mtxtal[1][10]={"YT's = $B"};
	static char mtxtam[1][58]={"Original ragged size specifications for dimenaion $M:$N$B"};
	static char mtxtan[1][49]={"Start of XT specifications for each dimension:$B"};
	static char mtxtao[1][53]={"Start of specifications for each ragged dimension:$B"};
	static char mtxtap[1][31]={"Final Ragged Specifications:$B"};
	static char mtxtaq[1][24]={"Dim. 1 to $M Indexes:$B"};
	static long mloc[12]={LTXTAA,LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,
	 LTXTAG,LTXTAH,LTXTAI,LTXTAJ,LTXTAK,LTXTAL};
	static long mactdg[3]={METDIG,6,MERET};
	static long mact[4]={MENTXT,0,METEXT,MERET};
	static long mactiv[4]={METEXT,MEIVEC,0,MERET};
	static long mactiv1[6]={MEMDA1,0,METEXT,MEIVEC,0,MERET};
	static long mactfv[4]={METEXT,MEFVEC,0,MERET};
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Eopt = &eopt[0] - 1;
	float *const Fdat = &fdat[0] - 1;
	long *const In = &in[0] - 1;
	long *const Intchk = &intchk[0] - 1;
	long *const Iopt = &iopt[0] - 1;
	long *const Knt = &knt[0] - 1;
	long *const Lr = &lr[0] - 1;
	long *const Lup = &lup[0] - 1;
	long *const Lx = &lx[0] - 1;
	long *const Mact = &mact[0] - 1;
	long *const Mactdg = &mactdg[0] - 1;
	long *const Mactfv = &mactfv[0] - 1;
	long *const Mactiv = &mactiv[0] - 1;
	long *const Mactiv1 = &mactiv1[0] - 1;
	long *const Mloc = &mloc[0] - 1;
	long *const Ndeg = &ndeg[0] - 1;
	long *const Ntab = &ntab[0] - 1;
	long *const Nup = &nup[0] - 1;
	float *const X = &x[0] - 1;
	float *const Xt = &xt[0] - 1;
	float *const Yt = &yt[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 2006, Math a la Carte, Inc.
	 *>> 2006-05-03 SILUPMD Krogh Started this debug routine
	 *
	 *--S replaces "?": ?ILUPMD, ?ILUPM, ?MESS
	 *
	 * Provide debugging information for Multidimensional Polynomial
	 * Interpolation with Look Up, silupm.
	 * Design/Code by Fred T. Krogh, Math a la Carte, Inc.
	 *
	 * *****************     Formal Arguments     ***************************
	 *
	 * These are exactly as for silupm.  See silupm.f for details.
	 *
	 * *****************     Variables Referenced     ***********************
	 *
	 * EOPT   Formal argument, see silupm.f.
	 * FDAT   Place to store floating point for error messages.
	 * I      Temporary index.
	 * IOPT   Formal argument, see silupm.f.
	 * IN     Array used to track current index in each dimension.
	 * INTCHK Array used for printing integers in error messages.
	 * IYT    Pointer to YT information for the current output.
	 * J      Index of a 0 in NTAB which is followed by user information
	 *        specifying the number of data points as a function of
	 *        previously used values from XT.
	 * K      Temporary index.
	 * KNT    Array used to keep count of number of values to output in
	 *        each dimension.
	 * LR     Arrray used to find ragged table information in each ragged
	 *        dimension.
	 * LTXTxx Parameter names of this form were generated by PMESS in making
	 *        up text for error messages and are used to locate various parts
	 *        of the message.
	 * LUP    Formal argument, see silupm.f.
	 * LX     Array giving start of XT table for each dimension.
	 * MACT   Array giving actions for printing error messages, see MESS.
	 * MAXDEG Parameter -- Gives maximum degree of polynomial interpolation.
	 * MECONT Parameter telling MESS an error message is to be continued.
	 * MEEMES Parameter telling MESS to print an error message.
	 * MEIVEC Parameter telling MESS to print an integer vector.
	 * MEFVEC Parameter telling MESS to print an floating point vector.
	 * MEMDA1 Parameter telling MESS that next item is integer data to be
	 *        made available for output.
	 * MERET  Parameter telling MESS that this ends the error message.
	 * MESS   Subroutine for printing error messages.
	 * METDIG Parameter for SMESS to set temporily the digits printed.
	 * NDEG   Formal argument, see silupm.f.
	 * NDIM   Formal argument, see silupm.f.
	 * NDIMI  Internal value for NDIM = number of dimensions.
	 * NTAB   Formal argument, see silupm.f.
	 * NTABM  = NTABXT + NDIMI -- NTAB(NTABM+I) is 1 + the index of the last
	 *        word in NTAB required for ragged table storage through dim. I.
	 * NTABXT = NDIMI+1 -- NTAB(NTABXT) = index of first ragged table, =1000
	 *        if there is no ragged table.  NTABN(NTABXT+I) = base address
	 *        accessing XT information.
	 * NUP    Array containing 0, unless an index change in this dimension
	 *        triggers a change in a ragged pointer.
	 * SMESS  Calls MESS and prints floating data in error messages.
	 * X      Formal argument, see silupm.f.
	 * XT     Formal argument, see silupm.f.
	 * Y      Formal argument, see silupm.f.
	 * YT     Formal argument, see silupm.f.
	 *
	 * *****************     Formal Variable Declarations     ***************
	 * */
 
	/* *****************    Parameters and Internal Varialbes ***************
	 *
	 *       MAXDEG and MAXDIM Should agree with values silupm.f */
 
	/* ************************ Error Message Stuff and Data ****************
	 *
	 * Parameter defined below are all defined in the error message program
	 * MESS and SMESS.
	 * */
 
 
	/* ********* Error message text ***************
	 * ********* Error message text ***************
	 *[Last 2 letters of Param. name]  [Text generating message.]
	 *AA Inputs to SILUPM, with NDIM = $I, EOPT dim. of $I.$E
	 *AB IOPT($I) = $I, is not a valid option.$E
	 *AC IOPT($I) = $I, requests an error estimate.$E
	 *AD IOPT($I) = $I, specifies$E
	 *AE IOPT($I) = $I, Store derivatives at $I, Max. total derivative $C
	 *   of $I,$E
	 *AF IOPT($I) = $I, requests Abs. & Rel. errors of: $F $F$E
	 *AG IOPT($I) = $I, request accuracy of $F$E
	 *AH IOPT($I) = $I, skips data points with YT(?) = $F$E
	 *AI LUP($I) = $I; it should be < 4.$E
	 *AJ Ragged table badly specified in dimension $I.$E
	 *AK Start of ragged table, NTAB($I) = $I; it must be 0.$E
	 *AL Middle of ragged table, NTAB($I) = $I; it must be >0.$E
	 *AM End of ragged table, NTAB($I) = $I; it must be -1.$E
	 *   $
	 *AN Extrpolation polynomial degrees:$B
	 *   $
	 *AO max derivatives per X(I):$B
	 *   $
	 *AP NDEG():$B
	 *   $
	 *AQ LUP():$B
	 *   $
	 *AR NTAB():$B
	 *   $
	 *AS X(): $B
	 *   $
	 *AT Dim. 1 XT's = $B
	 *   $
	 *AU Dim. 1 to $I XT's = $B
	 *   $
	 *AV Dim. $I XT's = $F + k * %F$E
	 *   $
	 *AW Dim. $I XT's = $B
	 *   $
	 *AX YT's = $B
	 *   $
	 *AY Original ragged size specifications for dimenaion $M:$N$B
	 *   $
	 *AZ Start of XT specifications for each dimension:$B
	 *   $
	 *BA Start of specifications for each ragged dimension:$B
	 *   $
	 *BB Final Ragged Specifications:$B
	 *   $
	 *BC Dim. 1 to $M Indexes:$B
	 *   $ */
	/* ********* End of code generated by PMESS ************** */
 
	/*                      1  2       3      4 */
	/*                        1       2  3      4 */
	/*                         1  2       3       4  5      5 */
 
	/* *****************     Start of executable code     *******************
	 * */
	ndimi = ndim;
	ntabxt = ndimi + 1;
	ntabm = ntabxt + ndimi;
	i = 0;
	Intchk[1] = ndimi;
	Intchk[2] = Iopt[2];
	Mact[2] = Mloc[1];
	smess( mactdg, (char*)mtxtaa,206, intchk, fdat );
L_10:
	;
	smess( mact, (char*)mtxtaa,206, intchk, fdat );
	if (i == 0)
	{
		if ((ndimi > MAXDIM) || (ndimi <= 0))
			{
			puts( "Code requires 0 < NDIM < 11." );
			exit(0);
			}
		if (Ntab[ntabxt] <= 0)
		{
			/* This block of code copied indented, else unchanged from silupm.f,
			 * except for code between lines starting with c### */
			Ntab[ntabxt] = 1000;
			Ntab[ntabxt + 1] = 1;
			for (i = 1; i <= ndimi; i++)
			{
				if (Lup[i] >= 4)
				{
					/*                       LUP(I) is out of range -- fatal error. */
					Intchk[1] = i;
					Intchk[2] = Lup[i];
					iiflg = -3;
				}
				l = Ntab[i];
				if (l > 0)
				{
					/*                             Table is not ragged up to this point. */
					if (i == ndimi)
						goto L_50;
					/*                       If table is ragged save pointer to ragged info. */
					if (Ntab[ndimi] < 0)
						Ntab[ntabm + i] = ntabm + ndimi;
					/*       Get pointer to start of XT data for next dimension */
					if (Lup[i] == 3)
					{
						Ntab[ntabxt + i + 1] = Ntab[ntabxt + i] +
						 2;
					}
					else
					{
						Ntab[ntabxt + i + 1] = Ntab[ntabxt + i] +
						 l;
					}
				}
				else if ((l == 0) || ((l != 1 - i) && ((i != ndimi) ||
				 (l <= -i))))
				{
					/*                                 Problem in specifying raggedness. */
					lt = ndimi;
					Intchk[1] = i;
					iiflg = -4;
					goto L_400;
				}
				else
				{
					j = Ntab[ntabm + i - 1];
					if (Ntab[i - 1] > 0)
					{
						Ntab[ntabxt] = -Ntab[i];
						/*       Get K = number of NTAB entries of extra data */
						k = Ntab[1];
						for (l = 2; l <= (-l); l++)
						{
							k *= Ntab[l];
						}
					}
					else
					{
						k = Ntab[j - 1];
					}
					if (Ntab[j] != 0)
					{
						lt = j;
						Intchk[1] = j;
						Intchk[2] = Ntab[j];
						iiflg = -5;
						goto L_400;
					}
					/*       Change data to be the partial sum of the original data. */
					lt = j + k;
 
 
					/*### This block of code added to that from silupm */
					Mactiv1[2] = i;
					Mactiv1[5] = lt - j + 1;
					mess( mactiv1, (char*)mtxtam,58, &Ntab[j] );
					/*### End of added block */
 
 
					for (l = j + 1; l <= lt; l++)
					{
						if (Ntab[l] <= 0)
						{
							/*                     Error, bad value inside ragged table info. */
							Intchk[1] = l;
							Intchk[2] = Ntab[l];
							iiflg = -6;
							goto L_400;
						}
						Ntab[l] += Ntab[l - 1];
					}
					if (i == ndimi)
					{
						lt += 1;
						if (Ntab[lt] == -1)
							goto L_50;
						/*                    Error -- No end tag where needed. */
						Intchk[1] = lt;
						Intchk[2] = Ntab[lt];
						iiflg = -7;
						goto L_400;
					}
					/*       Save index of 0th NTAB entry for extra data for next dim. */
					Ntab[ntabm + i] = j + k + 1;
					/*       Get pointer to start of XT data for next dimension */
					if (Lup[i] == 3)
					{
						Ntab[ntabxt + i + 1] = Ntab[ntabxt + i] +
						 2*k;
					}
					else
					{
						Ntab[ntabxt + i + 1] = Ntab[ntabxt + i] +
						 Ntab[j + k];
					}
				}
			}
		}
L_50:
		;
		/* End of code included from SILUPM */
		i = 2;
 
	}
	else if ((Intchk[2] <= 0) || (Intchk[2] > 6))
	{
		/*               A bad option value */
		exit(0);
	}
	else if (Intchk[2] == 2)
	{
		/*               Print out desired polynomial degrees for extrpolation. */
		Mactiv[3] = ndimi;
		mess( mactiv, (char*)mtxtab,35, &Iopt[i + 1] );
		i += ndimi;
	}
	else if (Intchk[2] == 3)
	{
		/*                Print out the derivative information */
		Mactiv[3] = ndimi;
		mess( mactiv, (char*)mtxtac,28, &Iopt[i + 3] );
		i += ndimi + 2;
	}
L_100:
	;
	i += 1;
	Intchk[1] = i;
	Intchk[2] = Iopt[i];
	Mact[2] = Mloc[Intchk[2] + 2];
	if (Intchk[2] == 0)
		goto L_200;
	/*             1   2    3    4    5    6 */
	switch (Intchk[2])
	{
		case 1: goto L_10;
		case 2: goto L_10;
		case 3: goto L_130;
		case 4: goto L_140;
		case 5: goto L_140;
		case 6: goto L_140;
	}
	Mact[2] = Mloc[2];
	goto L_10;
L_130:
	Intchk[3] = Iopt[i + 1];
	Intchk[4] = Iopt[i + 2];
	goto L_10;
L_140:
	Fdat[1] = Eopt[Iopt[i + 1]];
	if (Intchk[2] == 4)
		Fdat[2] = Eopt[Iopt[i + 1] + 1];
	i += 1;
	goto L_10;
 
L_200:
	;
	Mactiv[3] = ndimi;
	mess( mactiv, (char*)mtxtad,10, ndeg );
	mess( mactiv, (char*)mtxtae,9, lup );
	Mactiv[3] = ndimi + 1;
	mess( mactiv, (char*)mtxtaf,10, ntab );
	Mactiv[3] = ndimi;
	if (Ntab[ntabxt] == 0)
		goto L_300;
	mess( mactiv, (char*)mtxtan,49, &Ntab[ntabxt + 1] );
	if (Ntab[ntabxt] < 100)
	{
		mess( mactiv, (char*)mtxtao,53, &Ntab[ntabxt + ndimi + 1] );
		Mactiv[3] = lt - ntabxt - 2*ndimi;
		mess( mactiv, (char*)mtxtap,31, &Ntab[ntabxt + 2*ndimi + 1] );
	}
	Mactfv[3] = ndimi;
	smess( mactfv, (char*)mtxtag,8, intchk, x );
	Mactfv[3] = Ntab[1];
	smess( mactfv, (char*)mtxtah,17, intchk, xt );
 
	/* Set up for output of XT and YT */
	for (i = 1; i <= ndimi; i++)
	{
		Nup[i] = 0;
		In[i] = 1;
		Lx[i] = Ntab[ntabxt + i];
		if (Ntab[i] < 0)
		{
			Lr[i] = Ntab[2*ndimi + i];
			Knt[i] = Ntab[Lr[i] + 1] - Ntab[Lr[i]];
		}
		else
		{
			Knt[i] = Ntab[i];
		}
		Fdat[i] = Xt[Lx[i]];
	}
	iyt = 1;
L_240:
	;
	/*  Print XT, YT data for the current selection. */
	i = ndimi;
	Mactiv1[2] = i - 1;
	Mactiv1[5] = i - 1;
	mess( mactiv1, (char*)mtxtaq,24, in );
	Mactfv[3] = i - 1;
	Intchk[1] = i - 1;
	smess( mactfv, (char*)mtxtai,23, intchk, fdat );
	Mactfv[3] = Knt[i];
	Intchk[1] = i;
	if (Lup[i] == 3)
	{
		smess( &Mact[3], (char*)mtxtaj,29, intchk, &Xt[Lx[i]] );
	}
	else
	{
		smess( mactfv, (char*)mtxtak,18, intchk, &Xt[Lx[i]] );
	}
	smess( mactfv, (char*)mtxtal,10, intchk, &Yt[iyt] );
	iyt += Mactfv[3];
L_250:
	;
	if (Ntab[i] < 0)
		Nup[-Ntab[i]] = i;
	Fdat[i] = Xt[Lx[i]];
	i -= 1;
	/*               Test if we are done. */
	if (i == 0)
		goto L_300;
	if (Nup[i] != 0)
	{
		k = Nup[i];
		if (Lup[k] == 3)
		{
			Lx[k] += 2;
		}
		else
		{
			Lx[k] += Knt[k];
		}
		Fdat[k] = Xt[Lx[k]];
		Lr[k] += 1;
		Knt[k] = Ntab[Lr[k] + 1] - Ntab[Lr[k]];
	}
	if (In[i] < Knt[i])
	{
		if (Lup[i] == 3)
		{
			Fdat[i] = Xt[Lx[i]] + In[i]*Xt[Lx[i] + 1];
		}
		else
		{
			Fdat[i] = Xt[Lx[i] + In[i]];
		}
		In[i] += 1;
		goto L_240;
	}
	In[i] = 1;
	goto L_250;
 
	/* Restore default for number of digits before exit */
L_300:
	;
	smess( mactdg, (char*)mtxtaa,206, intchk, fdat );
	printf("End of output\n");
	return;
 
 
 
	/* More error message processing for errors that would show in silupm.f */
L_400:
	;
	/*  Restore state in NTAB. */
	if (Ntab[ntabxt] < 0)
	{
		for (i = 3*ndimi + 2; i <= (lt - 2); i++)
		{
			Ntab[i] = Ntab[i + 1] - Ntab[i];
		}
	}
	Ntab[ntabxt] = 0;
	Mact[2] = Mloc[7 - iiflg];
	mess( mact, (char*)mtxtaa,206, &Intchk[1] );
	i = 2;
	goto L_100;
} /* end of function */