share/hompack/fortran/dcpose.f

   1       SUBROUTINE DCPOSE(NDIM,N,QR,ALPHA,PIVOT,IERR,Y,SUM)
   2 C
   3 C SUBROUTINE  DCPOSE  IS A MODIFICATION OF THE ALGOL PROCEDURE
   4 C DECOMPOSE  IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST
   5 C SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS,
   6 C NUMER. MATH. 7 (1965) 269-276.
   7 C
   8       INTEGER NDIM,N,PIVOT(1)
   9       DOUBLE PRECISION QR(NDIM,1),ALPHA(N)
  10       INTEGER IERR,I,J,JBAR,K,KP1,NP1
  11       DOUBLE PRECISION BETA,SIGMA,ALPHAK,QRKK,Y(1),SUM(1)
  12       DOUBLE PRECISION DDOT
  13       IERR=0
  14       NP1=N+1
  15       DO 20 J=1,NP1
  16         SUM(J)=DDOT(N,QR(1,J),1,QR(1,J),1)
  17 20    PIVOT(J)=J
  18       DO 500 K=1,N
  19         SIGMA=SUM(K)
  20         JBAR=K
  21         KP1=K+1
  22         DO 40 J=KP1,NP1
  23           IF (SIGMA .GE. SUM(J)) GO TO 40
  24           SIGMA=SUM(J)
  25           JBAR=J
  26 40      CONTINUE
  27         IF (JBAR .EQ. K) GO TO 70
  28         I=PIVOT(K)
  29         PIVOT(K)=PIVOT(JBAR)
  30         PIVOT(JBAR)=I
  31         SUM(JBAR)=SUM(K)
  32         SUM(K)=SIGMA
  33         DO 50 I=1,N
  34           SIGMA=QR(I,K)
  35           QR(I,K)=QR(I,JBAR)
  36           QR(I,JBAR)=SIGMA
  37 50      CONTINUE
  38 C   END OF COLUMN INTERCHANGE.
  39 70      SIGMA=DDOT(N-K+1,QR(K,K),1,QR(K,K),1)
  40         IF (SIGMA .NE. 0.0) GO TO 60
  41         IERR=1
  42         RETURN
  43 60      IF (K .EQ. N) GO TO 500
  44         QRKK=QR(K,K)
  45         ALPHAK=-SQRT(SIGMA)
  46         IF (QRKK .LT. 0.0) ALPHAK=-ALPHAK
  47         ALPHA(K)=ALPHAK
  48         BETA=1.0/(SIGMA-QRKK*ALPHAK)
  49         QR(K,K)=QRKK-ALPHAK
  50         DO 80 J=KP1,NP1
  51 80      Y(J)=BETA*DDOT(N-K+1,QR(K,K),1,QR(K,J),1)
  52         DO 100 J=KP1,NP1
  53           DO 90 I=K,N
  54             QR(I,J)=QR(I,J)-QR(I,K)*Y(J)
  55 90        CONTINUE
  56           SUM(J)=SUM(J)-QR(K,J)**2
  57 100     CONTINUE
  58 500   CONTINUE
  59       ALPHA(N)=QR(N,N)
  60       RETURN
  61       END