/*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
#include "fcrt.h"
#include "shtcc.h"
#include
/* 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 */