SUBROUTINE DHTGEN (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 DHTGEN Krogh Added external statement. C>> 1994-11-11 DHTGEN Krogh Declared all vars. C>> 1994-10-20 DHTGEN Krogh Changes to use M77CON C>> 1987-08-19 DHTGEN Lawson Initial code. c--D 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: DAXPY, DDOT, DNRM2 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 DDOT, DNRM2 double precision U(*), UPARAM, C(*), DDOT, DNRM2 double precision 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.0D0, ZERO = 0.0D0) 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 = 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 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) + DDOT(NTERMS, U(IUL1),IUINC, C(I3),ICE) if (SUM .ne. ZERO) then FAC = SUM*BINV C(I2) = C(I2) + FAC*UPARAM call DAXPY(NTERMS, FAC, U(IUL1),IUINC, C(I3),ICE) endif 120 continue return end