C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C     DGRAM    10/10/74
C
C     PURPOSE
C     COMPUTE AN ORTHONORMAL BASIS FOR THE COLUMNS OF A MATRIX X USING
C     HOUSEHOLDER TRANSFORMATIONS.
C
C     USAGE
C     CALL DGRAM(X,M,N,EPS,IR,W)
C
C     ARGUMENTS
C     X   - INPUT M BY N MATRIX STORED COLUMNWISE (STORAGE MODE 0).  ON
C           RETURN CONTAINS AN M BY IR MATRIX CONTAINING AN ORTHONORMAL
C           BASIS FOR THE COLUMNS OF X STORED COLUMNWISE.
C           REAL*8
C     M   - NUMBER OF ROWS IN X
C           INTEGER*4
C     N   - NUMBER OF COLUMNS IN X
C           INTEGER*4
C     EPS - INPUT TOLERANCE FOR COMPUTING THE RANK OF X.  A REASONABLE
C           VALUE IS 1.D-13.
C           REAL*8
C     IR  - COMPUTED RANK OF X.  NUMBER OF COLUMNS OF X ON RETURN.
C           INTEGER*4
C     W   - WORKVECTOR OF LENGTH N.
C           REAL*8
C
C
      SUBROUTINE DGRAM(X,M,N,EPS,IR,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 X(M,N),W(1)
C
C     SCALE X
C
      DO 20 J=1,N
      SUM=0.D0
      DO 10 I=1,M
10    SUM=SUM+X(I,J)**2
20    W(J)=DSQRT(SUM)
      SMAX=0.D0
      DO 30 J=1,N
      IF(W(J).LT.SMAX) GO TO 30
      SMAX=W(J)
      JMAX=J
30    CONTINUE
      TOL=SMAX*EPS
      DO 40 J=1,N
      S=0.D0
      IF(W(J).GE.TOL) S=1.D0/W(J)
40    W(J)=S
      W(JMAX)=(1.D0+1.D-13)/SMAX
      DO 50 J=1,N
      S=W(J)
      DO 50 I=1,M
50    X(I,J)=X(I,J)*S
C
C     FACTOR TO LOWER TRIANGLE AND ORTHOGONAL MATRIX
C
      N1=N-1
      IR=N
      DO 150 K=1,N
C     CHOOSE COL FROM WHICH TO COMPUTE HOUSEHOLDER MATRIX
      SMAX=0.D0
      IMAX=K
      DO 102 J=K,N
C     COMPUTE THE LENGTH OF A COL FROM DIAGONAL DOWN
      SUM=0.D0
      DO 101 I=K,M
101   SUM=SUM+X(I,J)**2
      IF(SUM.LT.SMAX) GO TO 102
      SMAX=SUM
      IMAX=J
102   CONTINUE
C     PEMUTE COLS
      DO 103 I=1,M
      A=X(I,K)
      X(I,K)=X(I,IMAX)
103   X(I,IMAX)=A
      IF(DSQRT(SMAX).GT.EPS) GO TO 109
      IR=K-1
      DO 104 J=K,N
      DO 104 I=1,M
104   X(I,J)=0.D0
      GO TO 151
109   SUM=SMAX
      TEMP=SUM
      SUM=DSQRT(SUM)
      AKK=X(K,K)
      IF(AKK.LT.0.D0) SUM=-SUM
C     MAKE COL K INTO TRANSFORMATION VECTOR U
110   X(K,K)=AKK+SUM
C     SAVE DIVISOR OF THIS TRANSFORMATION
      CONS=TEMP+SUM*AKK
      W(K)=CONS
C     TRANSFORM THE REMAINING COLMNS
      IF(N1.LT.K) GO TO 150
120   DO 140 J=K,N1
C     PROJECT COL K ON COL J+1
      SUM=0.D0
      DO 130 I=K,M
130   SUM=SUM+X(I,K)*X(I,J+1)
      TEMP=SUM/CONS
C     SUBTRACT APPROPRIATE MULTIPLE OF COL K FROM COL J+1
      DO 140 I=K,M
140   X(I,J+1)=X(I,J+1)-TEMP*X(I,K)
150   CONTINUE
151   CONTINUE
C
C     THE TRIANGLE MINUS ITS DIAGONAL IS IN THE UPPER PART OF X
C     IT IS DISCARDED
C     CONVERT IMPLICIT FORM BACK TO EXPLICIT ORTHONORMAL BASIS FOR X
C     APPLY ORTHOGONAL TRANSFORMATIONS IN REVERSE ORDER TO UNIT VECTORS
C     IN REVERSE ORDER
C
      NCOLX=IR
      DO 210 K=1,NCOLX
      KK=NCOLX-K+1
C     APPLY TRANSFORMATION KK TO UNIT VECTOR KK
      TEMP=X(KK,KK)/W(KK)
      X(KK,KK)=1.D0-TEMP*X(KK,KK)
      IF(M.EQ.KK) GO TO 170
      KK1=KK+1
      DO 160 I=KK1,M
160   X(I,KK)=-TEMP*X(I,KK)
C     APPLY THE REMAINING TRANSFORMATION
170   KM=KK-1
      IF(KM.EQ.0) GO TO 220
      DO 180 I=1,KM
180   X(I,KK)=0.D0
      DO 210 L=1,KM
      LL=KK-L
      SUM=0.D0
      DO 190 I=LL,M
190   SUM=SUM+X(I,LL)*X(I,KK)
      TEMP=SUM/W(LL)
      DO 200 I=LL,M
200   X(I,KK)=X(I,KK)-TEMP*X(I,LL)
210   CONTINUE
220   RETURN
      END