/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:13 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dsvdrs s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_dsvdrs.h"
/*     program DRDSVDRS
 *>> 1996-06-25 DRDSVDRS   Krogh Changes for conversion to C.
 *>> 1996-05-28 DRDSVDRS   Krogh Added external statement.
 *>> 1994-10-19 DRDSVDRS   Krogh  Changes to use M77CON
 *>> 1992-03-18 DRDSVDRS   CLL Added "c" to "program" line above.
 *>> 1992-03-12 DRDSVDRS   CLL
 *     Demo driver for DSVDRS and DCOV3, Singular Value Decomposition
 *     and Covariance matrix.
 *
 *     The sample data was computed as
 *     y = 0.5 + 0.25 * sin(2*pi*x) + 0.125 * exp(-x)
 *     rounded to four decimal places.
 *     ------------------------------------------------------------------
 *--D replaces "?": DR?SVDRS, ?SVDRS, ?COV3, ?DOT, ?NRM2
 *     ------------------------------------------------------------------ */
		/* PARAMETER translations */
#define	MMAX	11
#define	NMAX	3
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, ierr, j, m, n;
	double a[NMAX][MMAX], b[MMAX], c[NMAX], dof, pi, rnorm, sing[NMAX],
	 stddev, var, work[2*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};
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const B = &b[0] - 1;
	double *const C = &c[0] - 1;
	double *const Sing = &sing[0] - 1;
	double *const Work = &work[0] - 1;
	double *const X = &x[0] - 1;
	double *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
 
 
	/*     ------------------------------------------------------------------ */
	m = MMAX;
	n = NMAX;
	printf(" DRDSVDRS..    Demo driver for DSVDRS and DCOV3.\n");
 
	pi = 4.0e0*atan( 1.0e0 );
	for (i = 1; i <= m; i++)
	{
		a[0][i - 1] = 1.0e0;
		a[1][i - 1] = sin( 2.0e0*pi*X[i] );
		a[2][i - 1] = exp( -X[i] );
		B[i] = Y[i];
	}
 
	dsvdrs( (double*)a, MMAX, m, n, b, MMAX, 1, sing, work );
 
	printf("\n Computed singular values:\n ");
	for (i = 1; i <= n; i++)
	{
		printf("%14.6g", Sing[i]);
	}
	printf("\n");
	printf("\n V matrix:\n");
	for (i = 1; i <= n; i++)
	{
		printf(" ");
		for (j = 1; j <= n; j++)
		{
			printf("%12.6f", a[j - 1][i - 1]);
		}
		printf("\n");
	}
 
	rnorm = dnrm2( m - n, &B[n + 1], 1 );
	dof = m - n;
	var = SQ(rnorm)/dof;
	stddev = sqrt( var );
	printf("\n Estimated Std. Dev. of data error =%10.6f\n", stddev);
 
	/*     Writing the singular value decomposition of A as
	 *              A = U * S * Transpose(V),
	 *     the solution of the least-squares problem A*c ~ b will be
	 *              c  = V * Inverse(S) *Transpose(U) * b.
	 *     We now have V in A(), S in SING(), and the product, Transpose(U)*b
	 *     in B().  So we can now compute c, as long as all singular values
	 *     are nonzero.
	 *     SING() contains the singular values as nonnegative
	 *     numbers in nonincreasing order, so we only need to test SING(N).
	 * */
	if (Sing[n] != 0.0e0)
	{
		for (i = 1; i <= n; i++)
		{
			B[i] /= Sing[i];
		}
		for (i = 1; i <= n; i++)
		{
			C[i] = ddot( n, &a[0][i - 1], MMAX, b, 1 );
		}
 
		printf(" Solution coefficients   =");
		for (j = 1; j <= n; j++)
		{
			printf("%10.6f", C[j]);
		}
		printf("\n");
	}
 
	/*                      Compute covariance matrix.
	 * */
	dcov3( (double*)a, MMAX, n, sing, var, work, &ierr );
	printf("\n After DCOV3, IERR = %4ld\n", ierr);
	if (ierr == 0)
	{
		printf("\n Covariance matrix:\n");
		for (i = 1; i <= n; i++)
		{
			printf(" ");
			for (j = 1; j <= n; j++)
			{
				printf("%14.6g", a[j - 1][i - 1]);
			}
			printf("\n");
		}
	}
	exit(0);
} /* end of function */