/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:08 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dblas2 s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_dblas2.h"
/*     program DRDBLAS2
 *>> 1996-06-18 DRDBLAS2 Krogh  Minor format change for C conversion.
 *>> 1994-10-19 DRDBLAS2 Krogh  Changes to use M77CON
 *>> 1991-11-27 DRDBLAS2 CLL
 *>> 1987-12-09 Lawson  Initial Code.
 *
 *     Demonstrates the use of BLAS subroutines DROTG, DROT, DAXPY,
 *     and DCOPY to implement an algorithm for solving a linear
 *     least squares problem using sequential accumulation of the
 *     data and Givens orthogonal transformations.
 *     YTAB() contains rounded values of -2 + 2*X + 3*Exp(-X)
 *     -----------------------------------------------------------------
 *--D replaces "?": DR?BLAS2, ?ROTG, ?ROT, ?AXPY, ?COPY
 *     ----------------------------------------------------------------- */
		/* PARAMETER translations */
#define	MC	3
#define	MC1	(MC + 1)
#define	MXY	11
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int ixy, j, nc1, _i, _r;
	double c, coef[MC], div, estsd, rg[MC1][MC1], s, w[MC1], x, y;
	static double zero[1];
	static double xtab[MXY]={0.0e0,.1e0,.2e0,.3e0,.4e0,.5e0,.6e0,.7e0,
	 .8e0,.9e0,1.0e0};
	static double ytab[MXY]={1.00e0,.91e0,.86e0,.82e0,.81e0,.82e0,
	 .85e0,.89e0,.95e0,1.02e0,1.10e0};
	static long nxy = MXY;
	static long nc = MC;
	static int _aini = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Coef = &coef[0] - 1;
	double *const W = &w[0] - 1;
	double *const Xtab = &xtab[0] - 1;
	double *const Ytab = &ytab[0] - 1;
	double *const Zero = &zero[0] - 1;
		/* end of OFFSET VECTORS */
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		Zero[1] = 0.0e0;
		_aini = 0;
	}
 
 
	/*     ----------------------------------------------------------------- */
	nc1 = nc + 1;
	dcopy( MC1*MC1, zero, 0, (double*)rg, 1 );
	for (ixy = 1; ixy <= nxy; ixy++)
	{
		x = Xtab[ixy];
		y = Ytab[ixy];
		/*                                      Build new row of [A:B] in W(). */
		W[1] = 1.0e0;
		W[2] = x;
		W[3] = exp( -x );
		W[4] = y;
		/*                                      Process W() into [R:G]. */
		for (j = 1; j <= nc; j++)
		{
			drotg( &rg[j - 1][j - 1], &W[j], &c, &s );
			drot( nc1 - j, &rg[j][j - 1], MC1, &W[j + 1], 1, c, s );
		}
		drotg( &rg[nc1 - 1][nc1 - 1], &W[nc1], &c, &s );
	}
	/*                                      Begin: Solve triangular system. */
	dcopy( nc, &rg[nc1 - 1][0], 1, coef, 1 );
	for (j = nc; j >= 1; j--)
	{
		div = rg[j - 1][j - 1];
		if (div == 0.0e0)
		{
			printf("ERROR:ZERO DIVISOR AT J =%2ld\n", j);
			exit(0);
		}
		Coef[j] /= div;
		daxpy( j - 1, -Coef[j], &rg[j - 1][0], 1, coef, 1 );
	}
	/*                                        End: Solve triangular system.
	 * */
	printf(" Solution: COEF() = ");
	for (j = 1; j <= nc; j++)
	{
		printf("%8.3f", Coef[j]);
	}
	printf("\n");
	estsd = fabs( rg[nc1 - 1][nc1 - 1] )/sqrt( (double)( nxy - nc ) );
	printf("\n Estimated Std. Dev. of data errors =%9.5f\n", estsd);
	exit(0);
} /* end of function */