SUBROUTINE SHTGEN (MODE,LPIVOT,L1,M,U,LDU,COLU,UPARAM, * C,LDC,NVC,COLC) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. c>> 1996-03-30 SHTGEN Krogh Added external statement. C>> 1994-11-11 SHTGEN Krogh Declared all vars. C>> 1994-10-20 SHTGEN Krogh Changes to use M77CON C>> 1987-08-19 SHTGEN Lawson Initial code. c--S replaces "?": ?HTGEN, ?AXPY, ?DOT, ?NRM2 c C Construction and/or application of a single Householder C transformation.. Q = I + U*(U**T)/b c where I is the MxM identity matrix, b is a scalar, and U is an c M-dimensional Householder vector. c All vectors are M-vectors but only the components in positions c LPIVOT, and L1 through M, will be referenced. c This version, identified by GEN at the end of its name, c has the GENerality to handle the options of the U-vector being c stored either as a column or row of a matrix, and the vectors in c C() may be either column or row vectors. c ------------------------------------------------------------------ C Subroutine arguments c C MODE [in] = 1 OR 2 When MODE = 1 this subr determines the c parameters for a Householder transformation and applies c the transformation to NVC vectors. When MODE = 2 this c subr applies a previously determined Householder c transformation. C LPIVOT [in] The index of the pivot element. C L1,M [in] If L1 .le. M elements LPIVOT and L1 through M will c be referenced. If L1 .gt. M the subroutine returns c immediately. This may be regarded c as performing an identity transformation. C U() [inout] Contains the "pivot" vector. Typically U() will be c a two-dimensional array in the calling program and the pivot c vector may be either a column or row in this array. c When MODE = 1 this is the vector from which Householder c parameters are to be determined. c When MODE = 2 this is the result from previous computation c with MODE = 1. c LDU [in] Leading dimensioning parameter for U() in the calling c program where U() is a two-dimensional array. Gives c storage spacing between elements in a row of U() when U() is c regarded as a two-dimensional array. C COLU [in] True means the pivot vector is a column of the 2-dim c array U(). Thus the successive elements of the pivot vector c are at unit storage spacing. c False means the pivot vector is a row of the 2-dim array U() c Thus the storage spacing between successive elements is LDU. c UPARAM [inout] Holds a value that supplements the contents c of U() to complete the definition of a c Householder transformation. Computed when MODE = 1 and c reused when MODE = 2. c UPARAM is the pivot component of the Householder U-vector. C C() [inout] On entry contains a set of NVC M-vectors to which a c Householder transformation is to be applied. c On exit contains the set of transformed vectors. c Typically in the calling program C() will be a 2-dim array c with leading dimensioning parameter LDC. C These vectors are the columns of an M x NVC matrix in C(,) if c COLC = true, and are rows of an NVC x M matrix in C(,) if c COLC = false. C LDC [in] Leading dimension of C(,). Require LDC .ge. M if c COLC = true. Require LDC .ge. NVC if COLC = false. C NVC [in] Number of vectors in C(,) to be transformed. c If NVC .le. 0 no reference will be made to the array C(,). c COLC [in] True means the transformations are to be applied to c columns of the array C(,). False means the transformations c are to be applied to rows of the array C(,). c ------------------------------------------------------------------ c Subprograms referenced: SAXPY, SDOT, SNRM2 c ------------------------------------------------------------------ c This code was originally developed by Charles L. Lawson and c Richard J. Hanson at Jet Propulsion Laboratory in 1973. The c original code was described and listed in the book, c c Solving Least Squares Problems c C. L. Lawson and R. J. Hanson c Prentice-Hall, 1974 c c Feb, 1985, C. L. Lawson & S. Y. Chan, JPL. Adapted code from the c Lawson & Hanson book to Fortran 77 for use in the JPL MATH77 c library. c Prefixing subprogram names with S or D for s.p. or d.p. versions. c Using generic names for intrinsic functions. c Adding calls to BLAS and MATH77 error processing subrs in some c program units. c July, 1987. CLL. Changed user interface so method of specifying c column/row storage options is more language-independent. C ------------------------------------------------------------------ external SDOT, SNRM2 real U(*), UPARAM, C(*), SDOT, SNRM2 real B, FAC, HOLD, VNORM, ONE, ZERO, SUM, BINV integer MODE, LPIVOT, L1, M, LDU, LDC, NVC integer IUPIV, IUL1, IUINC, IUL0 integer ICE, ICV, I2, I3, INCR, NTERMS, J logical COLU, COLC parameter (ONE = 1.0E0, ZERO = 0.0E0) C ------------------------------------------------------------------ if (0.ge.LPIVOT .or. LPIVOT.ge.L1 .or. L1.gt.M) return if(COLU) then IUPIV = LPIVOT IUL1 = L1 IUINC = 1 else IUPIV = 1 + LDU * (LPIVOT-1) IUL1 = 1 + LDU * (L1-1) IUINC = LDU endif c if( MODE .eq. 1) then C ****** CONSTRUCT THE TRANSFORMATION. ****** IUL0 = IUL1 - IUINC if(IUL0 .eq. IUPIV) then VNORM = SNRM2(M-L1+2, U(IUL0), IUINC) else HOLD = U(IUL0) U(IUL0) = U(IUPIV) VNORM = SNRM2(M-L1+2, U(IUL0), IUINC) U(IUL0) = HOLD endif c if (U(IUPIV) .gt. ZERO) VNORM = -VNORM UPARAM = U(IUPIV)-VNORM U(IUPIV) = VNORM endif C ****** Apply the transformation I + U*(U**T)/B to C. **** C if (NVC.LE.0) return B = UPARAM * U(IUPIV) c Here B .le. 0. If B = 0., return. if (B .eq. ZERO) return BINV = ONE / B c I2 = 1 - ICV + ICE*(LPIVOT-1) c INCR = ICE * (L1-LPIVOT) if(COLC) then ICE = 1 ICV = LDC I2 = LPIVOT - LDC INCR = L1 - LPIVOT else ICE = LDC ICV = 1 I2 = ICE*(LPIVOT-1) INCR = ICE*(L1-LPIVOT) endif c NTERMS = M-L1+1 do 120 J = 1,NVC I2 = I2 + ICV I3 = I2 + INCR SUM = UPARAM * C(I2) + SDOT(NTERMS, U(IUL1),IUINC, C(I3),ICE) if (SUM .ne. ZERO) then FAC = SUM*BINV C(I2) = C(I2) + FAC*UPARAM call SAXPY(NTERMS, FAC, U(IUL1),IUINC, C(I3),ICE) endif 120 continue return end