C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C     DH1      7/18/74
C
C     PURPOSE
C     FOR A GIVEN VECTOR W COMPUTE THE VECTOR U AND SCALARS B AND S SUCH
C     THAT THE VECTOR V=QW SATISFIES:
C     1) COMPONANTS 1 THROUGH L OF V ARE THE SAME AS W,
C     2) COMPONANT L+1 OF V IS S,
C     3) COMPONANTS L+2 THROUGH L+J+1 OF V ARE THE SAME AS W,
C     4) COMPONANTS L+J+2 THROUGH M OF V ARE ZEROES,
C     WHERE Q IS THE HOUSEHOLDER TRANSFORMATION MATRIX:
C     Q=I           IF B=0
C     Q=I+(1/B)UU'  OTHERWISE
C
C     USAGE
C     CALL DH1(L,J,M,W,B,S)
C
C     ARGUMENTS
C     L - EITHER POSITIVE OR ZERO
C         INTEGER*4
C     J - EITHER POSITIVE OR ZERO
C         INTEGER*4
C     M   LENGTH OF W.  L, J, AND M MUST SATISFY L+J<M.
C         INTEGER*4
C     W - INPUT VECTOR OF LENGTH M. CONTAINS U ON RETURN.
C         REAL*8
C     B - SCALAR DESCRIBED ABOVE.
C         REAL*8
C     S - SCALAR DESCRIBED ABOVE.
C         REAL*8
C
C     REFERENCE
C     HANSON,R.G. AND LAWSON,C.L. EXTENSIONS AND APPLICATIONS OF THE
C     HOUSEHOLDER ALGORITHM FOR SOLVING LINEAR LEAST SQUARES PROBLEMS.
C     MATHEMATICS OF COMPUTATION,23. (1969).
C
C
C
      SUBROUTINE DH1(L,J,M,W,B,S)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 W(1)
      DATA TOL/1.D-38/
      WLP1=W(L+1)
      C=WLP1**2
      I=L+J+2
20    IF(I.GT.M) GO TO 30
      C=C+W(I)**2
      I=I+1
      GO TO 20
30    CONTINUE
      IF(C.LE.TOL) C=0.D0
      S=-DSQRT(C)
      IF(WLP1.LT.0.D0) S=-S
      UP=WLP1-S
      IF(C.EQ.0.D0) UP=0.D0
      B=S*UP
      K=L+J+1
      DO 40 I=1,K
40    W(I)=0.D0
      W(L+1)=UP
      RETURN
      END
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DH2       9/14/74
C
C     PURPOSE
C     GIVEN A MATRIX C COMPUTE THE MATRIX V=QC FROM OUTPUT OF DH1.
C
C     USAGE
C     CALL DH2(L,J,M,W,B,S,C,N)
C
C     ARGUMENTS
C     L - AS FOR DH1
C     J - AS FOR DH1
C     M - AS FOR DH1
C     W - AS RETURNED BY DH1
C     S - AS RETURNED BY DH1
C     B - AS RETURNED BY DH1
C     C - INPUT M BY N MATRIX STORED COLUMNWISE (STORAGE MODE 0).
C         CONTAINS V ON RETURN.
C         REAL*8
C     N - NUMBER OF COLUMNS IN C
C         INTEGER*4
C
C
      SUBROUTINE DH2(L,J,M,W,B,S,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 W(1),C(M,N)
      IF(B.EQ.0.D0) RETURN
      DO 90 IC=1,N
      V=W(L+1)*C(L+1,IC)
      I=L+J+2
20    IF(I.GT.M) GO TO 30
      V=V+W(I)*C(I,IC)
      I=I+1
      GO TO 20
30    CONTINUE
      V=V/B
      C(L+1,IC)=C(L+1,IC)+W(L+1)*V
      I=L+J+2
80    IF(I.GT.M) GO TO 90
      C(I,IC)=C(I,IC)+W(I)*V
      I=I+1
      GO TO 80
90    CONTINUE
      RETURN
      END
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DH3      9/14/74
C
C     PURPOSE
C     GIVEN A MATRIX C COMPUTE THE MATRIX V=CQ FROM OUTPUT OF DH1.
C
C     USAGE
C     CALL DH3(L,J,M,W,B,S,C,N)
C
C     ARGUMENTS
C     L - AS FOR DH1
C     J - AS FOR DH1
C     M - AS FOR DH1
C     W - AS RETURNED BY DH1
C     B - AS RETURNED BY DH1
C     S - AS RETURNED BY DH1
C     C - INPUT N BY M MATRIX STORED COLUMNWISE (STORAGE MODE 0)
C         CONTAINS V ON RETURN
C         REAL*8
C     N - NUMBER OF ROWS IN C
C         INTEGER*4
C
C
      SUBROUTINE DH3(L,J,M,W,B,S,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 W(1),C(N,M)
      IF(B.EQ.0.D0) RETURN
      DO 90 IR=1,N
      V=W(L+1)*C(IR,L+1)
      I=L+J+2
20    IF(I.GT.M) GO TO 30
      V=V+W(I)*C(IR,I)
      I=I+1
      GO TO 20
30    CONTINUE
      V=V/B
      C(IR,L+1)=C(IR,L+1)+W(L+1)*V
      I=L+J+2
80    IF(I.GT.M) GO TO 90
      C(IR,I)=C(IR,I)+W(I)*V
      I=I+1
      GO TO 80
90    CONTINUE
      RETURN
      END