/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:55 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "shtcc.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ONE	1.0e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ shtcc(
long mode,
long lpivot,
long l1,
long m,
float u[],
float *uparam,
float c[],
long ldc,
long nvc)
{
	long int i, iul0, j, jcbase, jcpiv;
	float b, binv, fac, hold, sum, vnorm;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const C = &c[0] - 1;
	float *const U = &u[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1996-03-30 SHTCC  Krogh   Added external statement.
	 *>> 1994-11-11 SHTCC  Krogh   Declared all vars.
	 *>> 1994-10-20 SHTCC  Krogh  Changes to use M77CON
	 *>> 1987-08-19 SHTCC  Lawson  Initial code.
	 *--S replaces "?": ?HTCC, ?NRM2
	 *
	 *        Construction and/or application of a single Householder
	 *     transformation..     Q = I + U*(U**T)/b
	 *     where I is the MxM identity matrix, b is a scalar, and U is an
	 *     M-dimensional Householder vector.
	 *        All vectors are M-vectors but only the components in positions
	 *     LPIVOT, and L1 through M, will be referenced.
	 *        This version, identified by CC at the end of its name, is
	 *     specialized for the Column-Column case, i.e. U() is a vector or
	 *     a column of a matrix and C() is regarded a containing a set
	 *     of column vectors to which transformations will be applied.
	 *     ------------------------------------------------------------------
	 *                         Subroutine arguments
	 *
	 *     MODE  [in]  = 1 OR 2   When MODE = 1 this subr determines the
	 *           parameters for a Householder transformation and applies
	 *           the transformation to NVC vectors.  When MODE = 2 this
	 *           subr applies a previously determined Householder
	 *           transformation.
	 *     LPIVOT  [in]  The index of the pivot element.
	 *     L1,M  [in]  If L1 .le. M  elements LPIVOT and L1 through M will
	 *           be referenced.  If L1 .gt. M the subroutine returns
	 *           immediately.  This may be regarded
	 *           as performing an identity transformation.
	 *     U()  [inout]  Contains an M-dimensional vector with storage
	 *           spacing of 1 between elements.
	 *           When MODE = 1 this is the vector from which Householder
	 *           parameters are to be determined.
	 *           When MODE = 2 this is the result from previous computation
	 *           with MODE = 1.
	 *     UPARAM  [inout]  Holds a value that supplements the
	 *           contents of U() to complete the definition of a
	 *           Householder transformation.  Computed when MODE = 1 and
	 *           reused when MODE = 2.
	 *           UPARAM is the pivot component of the Householder U-vector.
	 *     C()  [inout]  On entry contains a set of NVC M-vectors to which a
	 *          Householder transformation is to be applied.
	 *          On exit contains the set of transformed vectors.
	 *          These vectors are the columns of an M x NVC matrix in C(,).
	 *     LDC  [in]  Leading dimension of C(,).  Require LDC .ge. M.
	 *     NVC  [in]  Number of vectors in C(,) to be transformed.
	 *           If NVC .le. 0 no reference will be made to the array C(,).
	 *     ------------------------------------------------------------------
	 *     Subprograms referenced: SNRM2
	 *     ------------------------------------------------------------------
	 *          This code was originally developed by Charles L. Lawson and
	 *     Richard J. Hanson at Jet Propulsion Laboratory in 1973.  The
	 *     original code was described and listed in the book,
	 *
	 *                  Solving Least Squares Problems
	 *                  C. L. Lawson and R. J. Hanson
	 *                  Prentice-Hall, 1974
	 *
	 *     Feb, 1985, C. L. Lawson & S. Y. Chan, JPL.  Adapted code from the
	 *     Lawson & Hanson book to Fortran 77 for use in the JPL MATH77
	 *     library.
	 *     Prefixing subprogram names with S or D for s.p. or d.p. versions.
	 *     Using generic names for intrinsic functions.
	 *     Adding calls to BLAS and MATH77 error processing subrs in some
	 *     program units.
	 *     July, 1987. CLL.  Changed user interface so method of specifying
	 *     column/row storage options is more language-independent.
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	if ((0 >= lpivot || lpivot >= l1) || l1 > m)
		return;
	if (mode == 1)
	{
		/*                            ****** CONSTRUCT THE TRANSFORMATION. ****** */
		iul0 = l1 - 1;
		if (iul0 == lpivot)
		{
			vnorm = snrm2( m - l1 + 2, &U[iul0], 1 );
		}
		else
		{
			hold = U[iul0];
			U[iul0] = U[lpivot];
			vnorm = snrm2( m - l1 + 2, &U[iul0], 1 );
			U[iul0] = hold;
		}
 
		if (U[lpivot] > ZERO)
			vnorm = -vnorm;
		*uparam = U[lpivot] - vnorm;
		U[lpivot] = vnorm;
	}
	/*            ****** Apply the transformation  I + U*(U**T)/B  to C. ****
	 * */
	if (nvc <= 0)
		return;
	b = *uparam*U[lpivot];
	/*                                 Here B .le. 0.  If B  =  0., return. */
	if (b == ZERO)
		return;
	binv = ONE/b;
	jcbase = 0;
 
	for (j = 1; j <= nvc; j++)
	{
		jcpiv = jcbase + lpivot;
		sum = C[jcpiv]**uparam;
		for (i = l1; i <= m; i++)
		{
			sum += C[jcbase + i]*U[i];
		}
		if (sum != ZERO)
		{
			fac = sum*binv;
			C[jcpiv] += fac**uparam;
			for (i = l1; i <= m; i++)
			{
				C[jcbase + i] += fac*U[i];
			}
		}
		jcbase += ldc;
	}
	return;
} /* end of function */