/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:14 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_sblas4 s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_sblas4.h"
/*     program DRSBLAS4
 *>> 1994-10-19 DRSBLAS4 Krogh  Changes to use M77CON
 *>> 1992-05-12 DRSBLAS4 CLL Improved consistency of various versions.
 *>> 1992-04-15 CLL
 *>> 1992-03-16 CLL
 *>> 1991-12-02 DRSBLAS4 CLL
 *>> 1991-07-31 DRSBLAS4 CLL
 *
 *     Demonstrates the use of BLAS subroutines SROTMG, SROTM, SAXPY,
 *     and SCOPY to implement an algorithm for solving a linear
 *     least squares problem using sequential accumulation of the
 *     data and Modified (Fast) Givens orthogonal transformations.
 *     YTAB() contains rounded values of -2 + 2*X + 3*Exp(-X)
 *
 *     ------------------------------------------------------------------
 *--S replaces "?": DR?BLAS4, ?AXPY, ?ROTMG, ?ROTM, ?COPY, ?VECP
 *     ------------------------------------------------------------------ */
		/* PARAMETER translations */
#define	MC	3
#define	MC1	(MC + 1)
#define	MC2	(MC + 2)
#define	MXY	11
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int ixy, j, nc1, nc2, _i, _r;
	float c[MC], div, estsd, param[5], rg[MC2][MC1], w[MC2], x, y;
	static float one[1], zero[1];
	static float xtab[MXY]={0.0e0,.1e0,.2e0,.3e0,.4e0,.5e0,.6e0,.7e0,
	 .8e0,.9e0,1.0e0};
	static float ytab[MXY]={1.00e0,.91e0,.86e0,.82e0,.81e0,.82e0,.85e0,
	 .89e0,.95e0,1.02e0,1.10e0};
	static long nxy = 11;
	static long nc = 3;
	static int _aini = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const C = &c[0] - 1;
	float *const One = &one[0] - 1;
	float *const Param = &param[0] - 1;
	float *const W = &w[0] - 1;
	float *const Xtab = &xtab[0] - 1;
	float *const Ytab = &ytab[0] - 1;
	float *const Zero = &zero[0] - 1;
		/* end of OFFSET VECTORS */
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		One[1] = 1.0e0;
		Zero[1] = 0.0e0;
		_aini = 0;
	}
 
 
	/*     ------------------------------------------------------------------ */
	printf(" DRSBLAS4..  Demo for SROTMG and SROTM.\n\n");
	nc1 = nc + 1;
	nc2 = nc1 + 1;
	scopy( MC1*MC2, zero, 0, (float*)rg, 1 );
	scopy( nc1, one, 0, &rg[nc2 - 1][0], 1 );
	for (ixy = 1; ixy <= nxy; ixy++)
	{
		x = Xtab[ixy];
		y = Ytab[ixy];
 
		/*          Set W() to be the next row of the least-squares problem.
		 * */
		W[1] = 1.0e0;
		W[2] = x;
		W[3] = expf( -x );
		W[4] = y;
		W[5] = 1.0e0;
 
		/*        Accumulate W() into the triangular matrix [R:G]
		 *        using Givens rotations.
		 * */
		for (j = 1; j <= nc1; j++)
		{
			srotmg( &rg[nc2 - 1][j - 1], &W[nc2], &rg[j - 1][j - 1],
			 W[j], param );
			if (j < nc1)
				srotm( nc1 - j, &rg[j][j - 1], MC1, &W[j + 1], 1,
				 param );
		}
 
	}
 
	/*               Solve the triangular system.
	 * */
	scopy( nc, &rg[nc1 - 1][0], 1, c, 1 );
	for (j = nc; j >= 1; j--)
	{
		div = rg[j - 1][j - 1];
		if (div == 0.0e0)
		{
			printf("ERROR:ZERO DIVISOR AT J =%ld \n", j);
			exit(0);
		}
		C[j] /= div;
		saxpy( j - 1, -C[j], &rg[j - 1][0], 1, c, 1 );
	}
 
	svecp( c, nc, "0 Solution: C()=" );
	printf(" \n");
	estsd = fabsf( rg[nc1 - 1][nc1 - 1] )*sqrtf( rg[nc2 - 1][nc1 - 1] )/
	 sqrtf( (float)( nxy - nc ) );
	printf("Estimated STD. DEV. of data errors =%g \n", estsd);
	exit(0);
} /* end of function */