/*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_dhfti s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_dhfti.h"
/*     program DRDHFTI
 *>> 2001-05-22 DRDHFTI  Krogh Minor change for making .f90 version.
 *>> 1996-07-03 DRDHFTI  Krogh  Special code for C conversion.
 *>> 1994-10-19 DRDHFTI  Krogh  Changes to use M77CON
 *>> 1987-12-09 DRDHFTI  Lawson  Initial Code.
 *                       Demo driver for DHFTI and DCOV2
 *--D replaces "?": DR?HFTI, ?HFTI, ?COV2
 *
 *     The sample data was computed as
 *     y = 0.5 + 0.25 * sin(2*pi*x) + 0.125 * exp(-x)
 *     rounded to four decimal places.
 *     ------------------------------------------------------------------
 *++ Code for .C. is active */
   long int k;
		/* PARAMETER translations */
#define	FOUR	4.0e0
#define	MMAX	11
#define	NMAX	3
#define	ONE	1.0e0
#define	TWO	2.0e0
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, ierr, ip[NMAX], j, krank;
	double a[NMAX][MMAX], c[MMAX], dof, pi, rnorm[1], stddev, var,
	 work[NMAX];
	static double x[MMAX]={0.0e0,0.1e0,0.2e0,0.3e0,0.4e0,0.5e0,0.6e0,
	 0.7e0,0.8e0,0.9e0,1.0e0};
	static double y[MMAX]={0.6250e0,0.7601e0,0.8401e0,0.8304e0,0.7307e0,
	 0.5758e0,0.4217e0,0.3243e0,0.3184e0,0.4039e0,0.5460e0};
	static long m = MMAX;
	static long n = NMAX;
	static long nc = 1;
	static double tau = 0.0e0;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const C = &c[0] - 1;
	long *const Ip = &ip[0] - 1;
	double *const Rnorm = &rnorm[0] - 1;
	double *const Work = &work[0] - 1;
	double *const X = &x[0] - 1;
	double *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
 
/*++ End */
	/*     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	pi = FOUR*atan( ONE );
	for (i = 1; i <= m; i++)
	{
		a[0][i - 1] = ONE;
		a[1][i - 1] = sin( TWO*pi*X[i] );
		a[2][i - 1] = exp( -X[i] );
		C[i] = Y[i];
	}
 
	dhfti( (double*)a, MMAX, m, n, c, MMAX, nc, tau, &krank, rnorm,
	 work, ip );
	dof = m - n;
	stddev = Rnorm[1]/sqrt( dof );
	var = SQ(stddev);
	printf(" Rank of linear system   =%4ld\n", krank);
	printf(" Std. Dev. of data error =%10.6f\n", stddev);
	printf(" Solution coefficients   =");
	for (j = 1; j <= n; j++)
	{
		printf("%10.6f", C[j]);
	}
	printf("\n");
 
	dcov2( (double*)a, MMAX, n, ip, var, &ierr );
	printf(" Error flag from DCOV2   =%4ld\n", ierr);
	printf(" Covariance matrix of computed coefficients:\n");
	printf(" \n");
	/*++ Code for ~.C. is inactive
	 *      do 30 I = 1,N
	 *         print '(1x,3(3x,2i3,g16.8))', (I,J,A(I,J),J=I,N)
	 *   30 continue
	 *++ Code for .C. is active */
   for (i = 0; i < n; i++){
      for (j = i; j < n; j+=3){
         for (k = j; k < (j < n - 3 ? j+3 : n); k++)
              printf("   %3ld%3ld%16.8g", i+1, k+1, a[k][i]);
         printf("\n");}
      }
	exit(0);
	/*++ End */
} /* end of function */