/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:10 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dilupm s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_dilupm.h"
		/* PARAMETER translations */
#define	H1	.08e0
#define	H2	.12e0
#define	ND1	50
#define	ND2	40
#define	NDEG1	8
#define	NDEG2	10
#define	NDIM	2
#define	NDIMT	(2*NDIM + 1)
#define	NEOPT	(2*(2*NDIM + NDEG1) + 1)
#define	NUMTAB	(ND1*ND2)
#define	NX	3
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, j, _i, _r;
	static long int ntab[NDIMT];
	double ans, eopt[NEOPT], x[NDIM], xt[2*NDIM], y;
	static double x1[NX]={2.7e0,.93e0,.765e0};
	static double x2[NX]={-.1e0,.05e0,.87e0};
	static long iopt[4]={0,NEOPT,1,0};
	static long ndeg[NDIM]={NDEG1,NDEG2};
	static long lup[NDIM]={3,3};
	static int _aini = 1;
	/* EQUIVALENCE translations */
	double _e0[2000];
	double *const yt = (double*)_e0;
	double (*const yt2)[ND2] = (double(*)[ND2])_e0;
	/* end of EQUIVALENCE translations */
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Eopt = &eopt[0] - 1;
	long *const Iopt = &iopt[0] - 1;
	long *const Lup = &lup[0] - 1;
	long *const Ndeg = &ndeg[0] - 1;
	long *const Ntab = &ntab[0] - 1;
	double *const X = &x[0] - 1;
	double *const X1 = &x1[0] - 1;
	double *const X2 = &x2[0] - 1;
	double *const Xt = &xt[0] - 1;
	double *const Yt = &yt[0] - 1;
		/* end of OFFSET VECTORS */
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		Ntab[1] = ND1;
		Ntab[2] = ND2;
		Ntab[3] = 0;
		_aini = 0;
	}
 
	/*>> 2001-05-22 DRDILM  Krogh Minor change for making .f90 version.
	 *>> 1994-10-19 DRDILM  Krogh  Changes to use M77CON
	 *>> 1992-03-04 DRDILM  Krogh  Initial Code.
	 *--D replaces "?": DR?ILM, ?ILUPM
	 * Demonstration driver for DILUPM.
	 * Given a table of sin(xy), for x = 0, .08, .16, ..., and y = 0, .12, .2
	 * ... interpolates for (x,y) = (2.7, -.1), (.93, .05), (.765, .87) and
	 * obtains an error estimate.
	 * */
	/*                       Set IOPT to get an error estimate. */
 
	/* Compute the XT and YT tables */
	Xt[1] = 0.e0;
	Xt[2] = H1;
	Xt[3] = 0.e0;
	Xt[4] = H2;
	for (i = 1; i <= ND1; i++)
	{
		for (j = 1; j <= ND2; j++)
		{
			yt2[i - 1][j - 1] = sin( (double)( (i - 1)*(j - 1) )*H1*
			 H2 );
		}
	}
	printf("IOP(1)    X1      X2          Y     Est. Error True Error\n");
	for (i = 1; i <= NX; i++)
	{
		X[1] = X1[i];
		X[2] = X2[i];
		ans = sin( X[1]*X[2] );
		dilupm( NDIM, x, &y, ntab, xt, yt, ndeg, lup, iopt, eopt );
		printf("%5ld%9.4f%9.4f%12.8f%10.2e%11.2e\n", Iopt[1], X1[i], X2[i], y, Eopt[1],
		   y - ans);
	}
	exit(0);
} /* end of function */