/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:52 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dhtgen.h" #include /* PARAMETER translations */ #define ONE 1.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dhtgen( long mode, long lpivot, long l1, long m, double u[], long ldu, LOGICAL32 colu, double *uparam, double c[], long ldc, long nvc, LOGICAL32 colc) { long int i2, i3, ice, icv, incr, iuinc, iul0, iul1, iupiv, j, nterms; double b, binv, fac, hold, sum, vnorm; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const C = &c[0] - 1; double *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 DHTGEN Krogh Added external statement. *>> 1994-11-11 DHTGEN Krogh Declared all vars. *>> 1994-10-20 DHTGEN Krogh Changes to use M77CON *>> 1987-08-19 DHTGEN Lawson Initial code. *--D replaces "?": ?HTGEN, ?AXPY, ?DOT, ?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 GEN at the end of its name, * has the GENerality to handle the options of the U-vector being * stored either as a column or row of a matrix, and the vectors in * C() may be either column or row vectors. * ------------------------------------------------------------------ * 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 the "pivot" vector. Typically U() will be * a two-dimensional array in the calling program and the pivot * vector may be either a column or row in this array. * 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. * LDU [in] Leading dimensioning parameter for U() in the calling * program where U() is a two-dimensional array. Gives * storage spacing between elements in a row of U() when U() is * regarded as a two-dimensional array. * COLU [in] True means the pivot vector is a column of the 2-dim * array U(). Thus the successive elements of the pivot vector * are at unit storage spacing. * False means the pivot vector is a row of the 2-dim array U() * Thus the storage spacing between successive elements is LDU. * 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. * Typically in the calling program C() will be a 2-dim array * with leading dimensioning parameter LDC. * These vectors are the columns of an M x NVC matrix in C(,) if * COLC = true, and are rows of an NVC x M matrix in C(,) if * COLC = false. * LDC [in] Leading dimension of C(,). Require LDC .ge. M if * COLC = true. Require LDC .ge. NVC if COLC = false. * NVC [in] Number of vectors in C(,) to be transformed. * If NVC .le. 0 no reference will be made to the array C(,). * COLC [in] True means the transformations are to be applied to * columns of the array C(,). False means the transformations * are to be applied to rows of the array C(,). * ------------------------------------------------------------------ * Subprograms referenced: DAXPY, DDOT, DNRM2 * ------------------------------------------------------------------ * 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 (colu) { iupiv = lpivot; iul1 = l1; iuinc = 1; } else { iupiv = 1 + ldu*(lpivot - 1); iul1 = 1 + ldu*(l1 - 1); iuinc = ldu; } if (mode == 1) { /* ****** CONSTRUCT THE TRANSFORMATION. ****** */ iul0 = iul1 - iuinc; if (iul0 == iupiv) { vnorm = dnrm2( m - l1 + 2, &U[iul0], iuinc ); } else { hold = U[iul0]; U[iul0] = U[iupiv]; vnorm = dnrm2( m - l1 + 2, &U[iul0], iuinc ); U[iul0] = hold; } if (U[iupiv] > ZERO) vnorm = -vnorm; *uparam = U[iupiv] - vnorm; U[iupiv] = vnorm; } /* ****** Apply the transformation I + U*(U**T)/B to C. **** * */ if (nvc <= 0) return; b = *uparam*U[iupiv]; /* Here B .le. 0. If B = 0., return. */ if (b == ZERO) return; binv = ONE/b; /* I2 = 1 - ICV + ICE*(LPIVOT-1) * INCR = ICE * (L1-LPIVOT) */ if (colc) { ice = 1; icv = ldc; i2 = lpivot - ldc; incr = l1 - lpivot; } else { ice = ldc; icv = 1; i2 = ice*(lpivot - 1); incr = ice*(l1 - lpivot); } nterms = m - l1 + 1; for (j = 1; j <= nvc; j++) { i2 += icv; i3 = i2 + incr; sum = *uparam*C[i2] + ddot( nterms, &U[iul1], iuinc, &C[i3], ice ); if (sum != ZERO) { fac = sum*binv; C[i2] += fac**uparam; daxpy( nterms, fac, &U[iul1], iuinc, &C[i3], ice ); } } return; } /* end of function */