SUBROUTINE DCPOSE(NDIM,N,QR,ALPHA,PIVOT,IERR,Y,SUM)
C
C SUBROUTINE  DCPOSE  IS A MODIFICATION OF THE ALGOL PROCEDURE
C DECOMPOSE  IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST
C SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS,
C NUMER. MATH. 7 (1965) 269-276.
C
      INTEGER NDIM,N,PIVOT(1)
      DOUBLE PRECISION QR(NDIM,1),ALPHA(N)
      INTEGER IERR,I,J,JBAR,K,KP1,NP1
      DOUBLE PRECISION BETA,SIGMA,ALPHAK,QRKK,Y(1),SUM(1)
      DOUBLE PRECISION DDOT
      IERR=0
      NP1=N+1
      DO 20 J=1,NP1
        SUM(J)=DDOT(N,QR(1,J),1,QR(1,J),1)
20    PIVOT(J)=J
      DO 500 K=1,N
        SIGMA=SUM(K)
        JBAR=K
        KP1=K+1
        DO 40 J=KP1,NP1
          IF (SIGMA .GE. SUM(J)) GO TO 40
          SIGMA=SUM(J)
          JBAR=J
40      CONTINUE
        IF (JBAR .EQ. K) GO TO 70
        I=PIVOT(K)
        PIVOT(K)=PIVOT(JBAR)
        PIVOT(JBAR)=I
        SUM(JBAR)=SUM(K)
        SUM(K)=SIGMA
        DO 50 I=1,N
          SIGMA=QR(I,K)
          QR(I,K)=QR(I,JBAR)
          QR(I,JBAR)=SIGMA
50      CONTINUE
C   END OF COLUMN INTERCHANGE.
70      SIGMA=DDOT(N-K+1,QR(K,K),1,QR(K,K),1)
        IF (SIGMA .NE. 0.0) GO TO 60
        IERR=1
        RETURN
60      IF (K .EQ. N) GO TO 500
        QRKK=QR(K,K)
        ALPHAK=-SQRT(SIGMA)
        IF (QRKK .LT. 0.0) ALPHAK=-ALPHAK
        ALPHA(K)=ALPHAK
        BETA=1.0/(SIGMA-QRKK*ALPHAK)
        QR(K,K)=QRKK-ALPHAK
        DO 80 J=KP1,NP1
80      Y(J)=BETA*DDOT(N-K+1,QR(K,K),1,QR(K,J),1)
        DO 100 J=KP1,NP1
          DO 90 I=K,N
            QR(I,J)=QR(I,J)-QR(I,K)*Y(J)
90        CONTINUE
          SUM(J)=SUM(J)-QR(K,J)**2
100     CONTINUE
500   CONTINUE
      ALPHA(N)=QR(N,N)
      RETURN
      END