[CONTACT]

[ABOUT]

[POLICY]

[ADVERTISE]

C ALGORITHM COLLECTED ALGORITHMS FROM

Found at: ftp.icm.edu.pl:70/packages/netlib/toms-2014-06-10/686

C      ALGORITHM 686, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 16, NO. 4, PP. 369-377.
                                                                        QRU00010
C DRIVER FOR UPDATING THE QR DECOMPOSITION OF A 4 BY 3 MATRIX.          QRU00020
C                                                                       QRU00030
      PARAMETER(LDA=4,LDR=3)                                            QRU00040
      INTEGER M,N,I,K,INFO,J,NMBIT                                      QRU00050
      INTEGER IPOS(0:5),IP(LDR)                                         QRU00060
      DOUBLE PRECISION A(LDA,LDR),Q(LDA,LDR),R(LDR,LDR),B(LDA,LDR)      QRU00070
      DOUBLE PRECISION WORK(2*LDA+2*LDR+1),W(200),V(LDA),U(LDR)         QRU00080
      DOUBLE PRECISION RCOND,DELTA                                      QRU00090
C                                                                       QRU00100
C     INITIALIZATION OF A,Q,R,M,N                                       QRU00110
C                                                                       QRU00120
      DATA Q/5D-1,5D-1,5D-1,5D-1,5D-1,5D-1,-5D-1,-5D-1,5D-1,-5D-1,      QRU00130
     %5D-1,-5D-1/                                                       QRU00140
      DATA R/1D0,0D0,0D0,-1D0,1D0,0D0,-1D0,-1D0,1D0/                    QRU00150
      DATA A/5D-1,5D-1,5D-1,5D-1,0D0,0D0,-1D0,-1D0,-5D-1,-1.5D0,        QRU00160
     %5D-1,-5D-1/                                                       QRU00170
      DATA M/4/,N/3/                                                    QRU00180
C                                                                       QRU00190
      WRITE(6,10)'A=Q*R, MATRIX A:'                                     QRU00200
10    FORMAT(/1X,70A)                                                   QRU00210
      CALL MATPRI(A,LDA,M,N)                                            QRU00220
      WRITE(6,*)'MATRIX Q:'                                             QRU00230
      CALL MATPRI(Q,LDA,M,N)                                            QRU00240
      WRITE(6,*)'MATRIX R:'                                             QRU00250
      CALL MATPRI(R,LDR,N,N)                                            QRU00260
C                                                                       QRU00270
C     DELETE COLUMN 2 OF A, UPDATE Q AND R                              QRU00280
C                                                                       QRU00290
      K=2                                                               QRU00300
      DO 20 I=K+1,N                                                     QRU00310
        DO 30 J=1,M                                                     QRU00320
          A(J,I-1)=A(J,I)                                               QRU00330
30      CONTINUE                                                        QRU00340
20    CONTINUE                                                          QRU00350
C                                                                       QRU00360
      INFO=1                                                            QRU00370
      WRITE(6,10)'DELETE COLUMN 2 OF A AND UPDATE Q AND R BY DDELC'     QRU00380
      CALL DDELC(Q,LDA,M,N,R,LDR,K,V,INFO)                              QRU00390
      WRITE(6,40)'DDELC',INFO                                           QRU00400
40    FORMAT(' ON RETURN FROM ',A5,' INFO=',I1)                         QRU00410
C                                                                       QRU00420
C     FOR ALL SUBROUTINES INFO=0 ON EXIT INDICATES SUCCESFUL TERMINATIONQRU00430
C                                                                       QRU00440
      WRITE(6,50)(V(I),I=1,M)                                           QRU00450
50    FORMAT(' COLUMN RECOMPUTED BY DDELC:',4F8.3)                      QRU00460
      N=N-1                                                             QRU00470
      WRITE(6,10)'UPDATED MATRIX Q:'                                    QRU00480
      CALL MATPRI(Q,LDA,M,N)                                            QRU00490
      WRITE(6,*)'UPDATED MATRIX R: '                                    QRU00500
      CALL MATPRI(R,LDR,N,N)                                            QRU00510
C                                                                       QRU00520
C     PRINT MAGNITUDE OF ELEMENT OF RESIDUAL MATRIX A-Q*R OF LARGEST    QRU00530
C     MAGNITUDE, AS WELL AS FOR MATRIX Q'*Q-I                           QRU00540
C                                                                       QRU00550
      CALL MAXRES(A,LDA,Q,R,LDR,M,N)                                    QRU00560
      CALL ORTCHK(Q,LDA,M,N,WORK)                                       QRU00570
C                                                                       QRU00580
C     DELETE ROW 3 OF A, UPDATE Q AND R                                 QRU00590
C                                                                       QRU00600
      K=3                                                               QRU00610
      DO 100 I=1,N                                                      QRU00620
        DO 110 J=K+1,M                                                  QRU00630
          A(J-1,I)=A(J,I)                                               QRU00640
110     CONTINUE                                                        QRU00650
        A(M,I)=0D0                                                      QRU00660
100   CONTINUE                                                          QRU00670
C                                                                       QRU00680
      WRITE(6,10)'DELETE ROW 3 OF A AND UPDATE Q AND R BY DDELR'        QRU00690
      CALL DDELR(Q,LDA,M,N,R,LDR,K,U,WORK,INFO)                         QRU00700
      WRITE(6,40)'DDELR',INFO                                           QRU00710
      WRITE(6,120)(U(I),I=1,N)                                          QRU00720
120   FORMAT(' ROW RECOMPUTED BY DDELR: ',4F8.3)                        QRU00730
      M=M-1                                                             QRU00740
      WRITE(6,10)'UPDATED MATRIX Q:'                                    QRU00750
      CALL MATPRI(Q,LDA,M,N)                                            QRU00760
      WRITE(6,*)'UPDATED MATRIX R:'                                     QRU00770
      CALL MATPRI(R,LDR,N,N)                                            QRU00780
C                                                                       QRU00790
C     PRINT MAGNITUDE OF ELEMENT OF RESIDUAL MATRIX A-Q*R OF LARGEST    QRU00800
C     MAGNITUDE, AS WELL AS FOR MATRIX Q'*Q-I                           QRU00810
C                                                                       QRU00820
      CALL MAXRES(A,LDA,Q,R,LDR,M,N)                                    QRU00830
      CALL ORTCHK(Q,LDA,M,N,WORK)                                       QRU00840
C                                                                       QRU00850
C     INSERT NEW FIRST COLUMN V INTO A, UPDATE Q AND R                  QRU00860
      K=1                                                               QRU00870
      DO 200 I=N,K,-1                                                   QRU00880
        DO 210 J=1,M                                                    QRU00890
          A(J,I+1)=A(J,I)                                               QRU00900
210     CONTINUE                                                        QRU00910
200   CONTINUE                                                          QRU00920
      DO 220 I=1,M                                                      QRU00930
        V(I)=I                                                          QRU00940
        A(I,K)=V(I)                                                     QRU00950
220   CONTINUE                                                          QRU00960
      N=N+1                                                             QRU00970
C                                                                       QRU00980
C     UPPER BOUND FOR RECIPROCAL CONDITION NUMBER                       QRU00990
C                                                                       QRU01000
      RCOND=1D-6                                                        QRU01010
C                                                                       QRU01020
      WRITE(6,230)(V(I),I=1,M)                                          QRU01030
230   FORMAT(/1X,'NEW 1ST COLUMN TO BE INSERTED INTO A BY DINSC:',4F8.3)QRU01040
      CALL DINSC(Q,LDA,M,N,R,LDR,K,V,RCOND,WORK,INFO)                   QRU01050
      WRITE(6,40)'DINSC',INFO                                           QRU01060
      WRITE(6,240)RCOND                                                 QRU01070
240   FORMAT(' COMPUTED RECIPROCAL CONDITION NUMBER BY DINSC:',E8.1)    QRU01080
      WRITE(6,10)'UPDATED MATRIX Q:'                                    QRU01090
      CALL MATPRI(Q,LDA,M,N)                                            QRU01100
      WRITE(6,*)'UPDATED MATRIX R:'                                     QRU01110
      CALL MATPRI(R,LDR,N,N)                                            QRU01120
C                                                                       QRU01130
C     PRINT MAGNITUDE OF ELEMENT OF RESIDUAL MATRIX A-Q*R OF LARGEST    QRU01140
C     MAGNITUDE, AS WELL AS FOR MATRIX Q'*Q-I                           QRU01150
C                                                                       QRU01160
      CALL MAXRES(A,LDA,Q,R,LDR,M,N)                                    QRU01170
      CALL ORTCHK(Q,LDA,M,N,WORK)                                       QRU01180
C                                                                       QRU01190
C     INSERT ROW BETWEEN ROWS 2 AND 3 OF A, UPDATE Q AND R              QRU01200
C                                                                       QRU01210
      K=3                                                               QRU01220
      DO 300 J=1,N                                                      QRU01230
        DO 310 I=M,K,-1                                                 QRU01240
          A(I+1,J)=A(I,J)                                               QRU01250
310     CONTINUE                                                        QRU01260
        U(J)=J                                                          QRU01270
        A(K,J)=U(J)                                                     QRU01280
300   CONTINUE                                                          QRU01290
C                                                                       QRU01300
      M=M+1                                                             QRU01310
      WRITE(6,320)(U(I),I=1,N)                                          QRU01320
320   FORMAT(/1X,'NEW 3RD ROW TO BE INSERTED INTO A BY DINSR:',4F8.3)   QRU01330
      CALL DINSR(Q,LDA,M,N,R,LDR,K,U,WORK,INFO)                         QRU01340
      WRITE(6,40)'DINSR',INFO                                           QRU01350
      WRITE(6,10)'UPDATED MATRIX Q:'                                    QRU01360
      CALL MATPRI(Q,LDA,M,N)                                            QRU01370
      WRITE(6,*)'UPDATED MATRIX R:'                                     QRU01380
      CALL MATPRI(R,LDR,N,N)                                            QRU01390
C                                                                       QRU01400
C     PRINT MAGNITUDE OF ELEMENT OF RESIDUAL MATRIX A-Q*R OF LARGEST    QRU01410
C     MAGNITUDE, AS WELL AS FOR MATRIX Q'*Q-I                           QRU01420
C                                                                       QRU01430
      CALL MAXRES(A,LDA,Q,R,LDR,M,N)                                    QRU01440
      CALL ORTCHK(Q,LDA,M,N,WORK)                                       QRU01450
C                                                                       QRU01460
C     ADD RANK ONE MATRIX U*V' TO A, UPDATE Q AND R                     QRU01470
C                                                                       QRU01480
      DO 400 I=1,M                                                      QRU01490
        V(I)=5D-1                                                       QRU01500
400   CONTINUE                                                          QRU01510
      V(3)=2D0                                                          QRU01520
      DO 410 I=1,N                                                      QRU01530
        U(I)=(-1)**(I+1)                                                QRU01540
410   CONTINUE                                                          QRU01550
C                                                                       QRU01560
      DO 420 J=1,N                                                      QRU01570
        DO 430 I=1,M                                                    QRU01580
          A(I,J)=A(I,J)+U(J)*V(I)                                       QRU01590
430     CONTINUE                                                        QRU01600
420   CONTINUE                                                          QRU01610
C                                                                       QRU01620
      WRITE(6,10)                                                       QRU01630
     %'RANK ONE MATRIX V*U'' ADDED TO A, Q AND R UPDATED BY DRNK1'      QRU01640
      WRITE(6,440)'V',(V(I),I=1,M)                                      QRU01650
440   FORMAT(/1X,'VECTOR ',1A,':',4F8.3)                                QRU01660
      WRITE(6,450)'U',(U(I),I=1,N)                                      QRU01670
450   FORMAT(' VECTOR ',1A,':',4F8.3)                                   QRU01680
      CALL DRNK1(Q,LDA,M,N,R,LDR,U,V,WORK,INFO)                         QRU01690
      WRITE(6,460)'DRNK1',INFO                                          QRU01700
C                                                                       QRU01710
C     INFO=1 ON EXIT INDICATES THAT V ON ENTRY WAS IN SPAN(A)           QRU01720
C                                                                       QRU01730
460   FORMAT(/1X,'ON RETURN FROM ',A5,' INFO=',I1)                      QRU01740
                                                                        QRU01750
      WRITE(6,10)'UPDATED MATRIX Q:'                                    QRU01760
      CALL MATPRI(Q,LDA,M,N)                                            QRU01770
      WRITE(6,*)'UPDATED MATRIX R:'                                     QRU01780
      CALL MATPRI(R,LDR,N,N)                                            QRU01790
C                                                                       QRU01800
C     PRINT MAGNITUDE OF ELEMENT OF RESIDUAL MATRIX A-Q*R OF LARGEST    QRU01810
C     MAGNITUDE, AS WELL AS FOR MATRIX Q'*Q-I. PRINT A.                 QRU01820
C                                                                       QRU01830
      CALL MAXRES(A,LDA,Q,R,LDR,M,N)                                    QRU01840
      CALL ORTCHK(Q,LDA,M,N,WORK)                                       QRU01850
C                                                                       QRU01860
C     RANK REVEALING QR FACTORIZATION, INITIALIZE PERMUTATION VECTOR    QRU01870
C                                                                       QRU01880
      DO 500 J=1,N                                                      QRU01890
        IP(J)=J                                                         QRU01900
500   CONTINUE                                                          QRU01910
      K=N                                                               QRU01920
      NMBIT=2                                                           QRU01930
      DELTA=10                                                          QRU01940
C                                                                       QRU01950
      WRITE(6,10)'RANK REVEALING QR FACTORIZATION BY DRRPM'             QRU01960
      CALL DRRPM(Q,LDA,M,N,R,LDR,K,IP,DELTA,NMBIT,IPOS,WORK,INFO)       QRU01970
      WRITE(6,40)'DRRPM',INFO                                           QRU01980
      WRITE(6,510)DELTA                                                 QRU01990
510   FORMAT(' DELTA ON EXIT FROM DRRPM:',E8.1)                         QRU02000
      WRITE(6,520)(IPOS(J),J=0,NMBIT)                                   QRU02010
520   FORMAT(' POSITION OF ELEMENTS OF MAX MAGNITUDE OF SUCCESSIVELY ', QRU02020
     %'COMPUTED',/,' SINGULAR VECTORS BY DRRPM:',5I3)                   QRU02030
      WRITE(6,530)(IP(J),J=1,N)                                         QRU02040
530   FORMAT(' DRRPM DETERMINED QR FACTORIZATION OF MATRIX ',           QRU02050
     %'A(*,IP(J)), WHERE',/1X,' IP(1), IP(2), ... = ',10I3)             QRU02060
      WRITE(6,10)'MATRIX Q FOR COLUMN PERMUTED MATRIX A:'               QRU02070
      CALL MATPRI(Q,LDA,M,N)                                            QRU02080
      WRITE(6,*)'MATRIX R FOR COLUMN PERMUTED MATRIX A:'                QRU02090
      CALL MATPRI(R,LDR,N,N)                                            QRU02100
C                                                                       QRU02110
C     PERMUTE COLUMNS OF A ACCORDING TO IP, STORE IN B                  QRU02120
C                                                                       QRU02130
      DO 540 J=1,N                                                      QRU02140
        DO 550 K=1,M                                                    QRU02150
          B(K,J)=A(K,IP(J))                                             QRU02160
550     CONTINUE                                                        QRU02170
540   CONTINUE                                                          QRU02180
C                                                                       QRU02190
C     PRINT MAGNITUDE OF ELEMENT OF RESIDUAL MATRIX B-Q*R OF LARGEST    QRU02200
C     MAGNITUDE, AS WELL AS FOR MATRIX Q'*Q-I                           QRU02210
C                                                                       QRU02220
      CALL MAXRES(B,LDA,Q,R,LDR,M,N)                                    QRU02230
      CALL ORTCHK(Q,LDA,M,N,WORK)                                       QRU02240
C                                                                       QRU02250
      WRITE(6,10)'MATRIX A AFTER COLUMN PERMUTATION:'                   QRU02260
      CALL MATPRI(B,LDA,M,N)                                            QRU02270
C                                                                       QRU02280
      STOP                                                              QRU02290
      END                                                               QRU02300
C---------------------------------------------------------------------
C SUBROUTINES FOR COMPUTING RESIDUALS AND OUTPUT
C---------------------------------------------------------------------
      SUBROUTINE MATPRI(A,LDA,M,N)
C
C     MATPRI PRINTS MATRIX A
C
      INTEGER LDA,M,N,IROW,ICOL
      DOUBLE PRECISION A(LDA,N)
C
      DO 10 IROW=1,M
        WRITE(6,20)(A(IROW,ICOL), ICOL=1,N)
10    CONTINUE
      WRITE(6,*)
      RETURN
20    FORMAT(4F8.3)
      END
C--------------------------------------------------------------------
      SUBROUTINE ORTCHK(Q,LDQ,M,N,WK)
C
C     ORTCHK COMPUTES Q'*Q-I FOR THE M BY N MATRIX Q.
C
      INTEGER LDQ,M,N,IROW,ICOL,K
      DOUBLE PRECISION Q(LDQ,N),SUM,MX
C
      MX=0D0
      DO 10 IROW=1,N
        DO 20 ICOL=1,N
          SUM=0D0
          DO 30 K=1,M
            SUM=SUM+Q(K,IROW)*Q(K,ICOL)
30        CONTINUE
          IF(IROW.EQ.ICOL)SUM=SUM-1D0
          IF(MX.LT.DABS(SUM))MX=DABS(SUM)
20      CONTINUE
10    CONTINUE
C
      WRITE(6,40)MX
40    FORMAT(' ABS( ELEMENT OF LARGEST MAGNITUDE OF Q''*Q-I ):',
     %E8.1,/1X)
      RETURN
      END
C--------------------------------------------------------------------
      SUBROUTINE MAXRES(A,LDA,Q,R,LDR,M,N)
C
      DOUBLE PRECISION A(LDA,N),Q(LDA,M),R(LDR,N),SUM,MX
      MX=0D0
      DO 10 I=1,M
        DO 20 K=1,N
          SUM=A(I,K)
          DO 30 J=1,K
            SUM=SUM-Q(I,J)*R(J,K)
30        CONTINUE
          IF(DABS(SUM).GT.MX)MX=DABS(SUM)
20      CONTINUE
10    CONTINUE
      WRITE(6,40)MX
      RETURN
40    FORMAT(' ABS( ELEMENT OF LARGEST MAGNITUDE OF A-Q*R ):',E9.1)
      END
*
C Subject: subroutines for driver for subroutine testing
*
C---------------------------------------------------------------------
C SUBPROGRAMS FROM LINPACK FOR DRIVER OF UPDATE.FORTRAN
C
C SUBROUTINES: DQRDC, DSVDC
C BLAS: DAXPYX, DDOTX, DNRM2X, DROTGX, DSCALX, DSWAPX
C
C TO EACH NAME OF A LINPACK BLAS A TRAILING X HAS BEEN ATTACHED IN
C ORDER TO DISTINGUISH THESE BLAS FROM THE BLAS FOR THE UPDATING
C SUBROUTINES.
C---------------------------------------------------------------------
      SUBROUTINE DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
      INTEGER LDX,N,P,LDU,LDV,JOB,INFO
      DOUBLE PRECISION X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1)
C
C
C     DSVDC IS A SUBROUTINE TO REDUCE A DOUBLE PRECISION NXP MATRIX X
C     BY ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE
C     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE
C     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
C     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
C
C     ON ENTRY
C
C         X         DOUBLE PRECISION(LDX,P), WHERE LDX.GE.N.
C                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
C                   DECOMPOSITION IS TO BE COMPUTED.  X IS
C                   DESTROYED BY DSVDC.
C
C         LDX       INTEGER.
C                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.
C
C         N         INTEGER.
C                   N IS THE NUMBER OF COLUMNS OF THE MATRIX X.
C
C         P         INTEGER.
C                   P IS THE NUMBER OF ROWS OF THE MATRIX X.
C
C         LDU       INTEGER.
C                   LDU IS THE LEADING DIMENSION OF THE ARRAY U.
C                   (SEE BELOW).
C
C         LDV       INTEGER.
C                   LDV IS THE LEADING DIMENSION OF THE ARRAY V.
C                   (SEE BELOW).
C
C         WORK      DOUBLE PRECISION(N).
C                   WORK IS A SCRATCH ARRAY.
C
C         JOB       INTEGER.
C                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR
C                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB
C                   WITH THE FOLLOWING MEANING
C
C                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR
C                                  VECTORS.
C                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS
C                                  IN U.
C                        A.GE.2    RETURN THE FIRST MIN(N,P) SINGULAR
C                                  VECTORS IN U.
C                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR
C                                  VECTORS.
C                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS
C                                  IN V.
C
C     ON RETURN
C
C         S         DOUBLE PRECISION(MM), WHERE MM=MIN(N+1,P).
C                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
C                   SINGULAR VALUES OF X ARRANGED IN DESCENDING
C                   ORDER OF MAGNITUDE.
C
C         E         DOUBLE PRECISION(P).
C                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE
C                   DISCUSSION OF INFO FOR EXCEPTIONS.
C
C         U         DOUBLE PRECISION(LDU,K), WHERE LDU.GE.N.  IF
C                                   JOBA.EQ.1 THEN K.EQ.N, IF JOBA.GE.2
C                                   THEN K.EQ.MIN(N,P).
C                   U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P
C                   OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
C                   IN THE SUBROUTINE CALL.
C
C         V         DOUBLE PRECISION(LDV,P), WHERE LDV.GE.P.
C                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C                   V IS NOT REFERENCED IF JOB.EQ.0.  IF P.LE.N,
C                   THEN V MAY BE IDENTIFIED WITH X IN THE
C                   SUBROUTINE CALL.
C
C         INFO      INTEGER.
C                   THE SINGULAR VALUES (AND THEIR CORRESPONDING
C                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
C                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF
C                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
C                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX
C                   B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
C                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
C                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
C                   IS THE TRANSPOSE OF U).  THUS THE SINGULAR
C                   VALUES OF X AND B ARE THE SAME.
C
C     LINPACK. THIS VERSION DATED 03/19/79 .
C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C     DSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
C
C     EXTERNAL DROT
C     BLAS DAXPYX,DDOTX,DSCALX,DSWAPX,DNRM2X,DROTGX
C     FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT
C
C     INTERNAL VARIABLES
C
      INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
     *        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
      DOUBLE PRECISION DDOTX,T,R
      DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,DNRM2X,SCALE,SHIFT,SL,SM,SN,
     *                 SMM1,T1,TEST,ZTEST
      LOGICAL WANTU,WANTV
C
C
C     SET THE MAXIMUM NUMBER OF ITERATIONS.
C
      MAXIT = 30
C
C     DETERMINE WHAT IS TO BE COMPUTED.
C
      WANTU = .FALSE.
      WANTV = .FALSE.
      JOBU = MOD(JOB,100)/10
      NCU = N
      IF (JOBU .GT. 1) NCU = MIN0(N,P)
      IF (JOBU .NE. 0) WANTU = .TRUE.
      IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
C
C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
C
      INFO = 0
      NCT = MIN0(N-1,P)
      NRT = MAX0(0,MIN0(P-2,N))
      LU = MAX0(NCT,NRT)
      IF (LU .LT. 1) GO TO 170
      DO 160 L = 1, LU
         LP1 = L + 1
         IF (L .GT. NCT) GO TO 20
C
C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
C           PLACE THE L-TH DIAGONAL IN S(L).
C
            S(L) = DNRM2X(N-L+1,X(L,L),1)
            IF (S(L) .EQ. 0.0D0) GO TO 10
               IF (X(L,L) .NE. 0.0D0) S(L) = DSIGN(S(L),X(L,L))
               CALL DSCALX(N-L+1,1.0D0/S(L),X(L,L),1)
               X(L,L) = 1.0D0 + X(L,L)
   10       CONTINUE
            S(L) = -S(L)
   20    CONTINUE
         IF (P .LT. LP1) GO TO 50
         DO 40 J = LP1, P
            IF (L .GT. NCT) GO TO 30
            IF (S(L) .EQ. 0.0D0) GO TO 30
C
C              APPLY THE TRANSFORMATION.
C
               T = -DDOTX(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
               CALL DAXPYX(N-L+1,T,X(L,L),1,X(L,J),1)
   30       CONTINUE
C
C           PLACE THE L-TH ROW OF X INTO  E FOR THE
C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
C
            E(J) = X(L,J)
   40    CONTINUE
   50    CONTINUE
         IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
C
C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
C           MULTIPLICATION.
C
            DO 60 I = L, N
               U(I,L) = X(I,L)
   60       CONTINUE
   70    CONTINUE
         IF (L .GT. NRT) GO TO 150
C
C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
C           L-TH SUPER-DIAGONAL IN E(L).
C
            E(L) = DNRM2X(P-L,E(LP1),1)
            IF (E(L) .EQ. 0.0D0) GO TO 80
               IF (E(LP1) .NE. 0.0D0) E(L) = DSIGN(E(L),E(LP1))
               CALL DSCALX(P-L,1.0D0/E(L),E(LP1),1)
               E(LP1) = 1.0D0 + E(LP1)
   80       CONTINUE
            E(L) = -E(L)
            IF (LP1 .GT. N .OR. E(L) .EQ. 0.0D0) GO TO 120
C
C              APPLY THE TRANSFORMATION.
C
               DO 90 I = LP1, N
                  WORK(I) = 0.0D0
   90          CONTINUE
               DO 100 J = LP1, P
                  CALL DAXPYX(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
  100          CONTINUE
               DO 110 J = LP1, P
                  CALL DAXPYX(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
  110          CONTINUE
  120       CONTINUE
            IF (.NOT.WANTV) GO TO 140
C
C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
C              BACK MULTIPLICATION.
C
               DO 130 I = LP1, P
                  V(I,L) = E(I)
  130          CONTINUE
  140       CONTINUE
  150    CONTINUE
  160 CONTINUE
  170 CONTINUE
C
C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
C
      M = MIN0(P,N+1)
      NCTP1 = NCT + 1
      NRTP1 = NRT + 1
      IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
      IF (N .LT. M) S(M) = 0.0D0
      IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
      E(M) = 0.0D0
C
C     IF REQUIRED, GENERATE U.
C
      IF (.NOT.WANTU) GO TO 300
         IF (NCU .LT. NCTP1) GO TO 200
         DO 190 J = NCTP1, NCU
            DO 180 I = 1, N
               U(I,J) = 0.0D0
  180       CONTINUE
            U(J,J) = 1.0D0
  190    CONTINUE
  200    CONTINUE
         IF (NCT .LT. 1) GO TO 290
         DO 280 LL = 1, NCT
            L = NCT - LL + 1
            IF (S(L) .EQ. 0.0D0) GO TO 250
               LP1 = L + 1
               IF (NCU .LT. LP1) GO TO 220
               DO 210 J = LP1, NCU
                  T = -DDOTX(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
                  CALL DAXPYX(N-L+1,T,U(L,L),1,U(L,J),1)
  210          CONTINUE
  220          CONTINUE
               CALL DSCALX(N-L+1,-1.0D0,U(L,L),1)
               U(L,L) = 1.0D0 + U(L,L)
               LM1 = L - 1
               IF (LM1 .LT. 1) GO TO 240
               DO 230 I = 1, LM1
                  U(I,L) = 0.0D0
  230          CONTINUE
  240          CONTINUE
            GO TO 270
  250       CONTINUE
               DO 260 I = 1, N
                  U(I,L) = 0.0D0
  260          CONTINUE
               U(L,L) = 1.0D0
  270       CONTINUE
  280    CONTINUE
  290    CONTINUE
  300 CONTINUE
C
C     IF IT IS REQUIRED, GENERATE V.
C
      IF (.NOT.WANTV) GO TO 350
         DO 340 LL = 1, P
            L = P - LL + 1
            LP1 = L + 1
            IF (L .GT. NRT) GO TO 320
            IF (E(L) .EQ. 0.0D0) GO TO 320
               DO 310 J = LP1, P
                  T = -DDOTX(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
                  CALL DAXPYX(P-L,T,V(LP1,L),1,V(LP1,J),1)
  310          CONTINUE
  320       CONTINUE
            DO 330 I = 1, P
               V(I,L) = 0.0D0
  330       CONTINUE
            V(L,L) = 1.0D0
  340    CONTINUE
  350 CONTINUE
C
C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
C
      MM = M
      ITER = 0
  360 CONTINUE
C
C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
C
C     ...EXIT
         IF (M .EQ. 0) GO TO 620
C
C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
C        FLAG AND RETURN.
C
         IF (ITER .LT. MAXIT) GO TO 370
            INFO = M
C     ......EXIT
            GO TO 620
  370    CONTINUE
C
C        THIS SECTION OF THE PROGRAM INSPECTS FOR
C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
C        COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
C
C           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
C           KASE = 2     IF S(L) IS NEGLIGIBLE AND L.LT.M
C           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
C                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
C           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
C
         DO 390 LL = 1, M
            L = M - LL
C        ...EXIT
            IF (L .EQ. 0) GO TO 400
            TEST = DABS(S(L)) + DABS(S(L+1))
            ZTEST = TEST + DABS(E(L))
            IF (ZTEST .NE. TEST) GO TO 380
               E(L) = 0.0D0
C        ......EXIT
               GO TO 400
  380       CONTINUE
  390    CONTINUE
  400    CONTINUE
         IF (L .NE. M - 1) GO TO 410
            KASE = 4
         GO TO 480
  410    CONTINUE
            LP1 = L + 1
            MP1 = M + 1
            DO 430 LLS = LP1, MP1
               LS = M - LLS + LP1
C           ...EXIT
               IF (LS .EQ. L) GO TO 440
               TEST = 0.0D0
               IF (LS .NE. M) TEST = TEST + DABS(E(LS))
               IF (LS .NE. L + 1) TEST = TEST + DABS(E(LS-1))
               ZTEST = TEST + DABS(S(LS))
               IF (ZTEST .NE. TEST) GO TO 420
                  S(LS) = 0.0D0
C           ......EXIT
                  GO TO 440
  420          CONTINUE
  430       CONTINUE
  440       CONTINUE
            IF (LS .NE. L) GO TO 450
               KASE = 3
            GO TO 470
  450       CONTINUE
            IF (LS .NE. M) GO TO 460
               KASE = 1
            GO TO 470
  460       CONTINUE
               KASE = 2
               L = LS
  470       CONTINUE
  480    CONTINUE
         L = L + 1
C
C        PERFORM THE TASK INDICATED BY KASE.
C
         GO TO (490,520,540,570), KASE
C
C        DEFLATE NEGLIGIBLE S(M).
C
  490    CONTINUE
            MM1 = M - 1
            F = E(M-1)
            E(M-1) = 0.0D0
            DO 510 KK = L, MM1
               K = MM1 - KK + L
               T1 = S(K)
               CALL DROTGX(T1,F,CS,SN)
               S(K) = T1
               IF (K .EQ. L) GO TO 500
                  F = -SN*E(K-1)
                  E(K-1) = CS*E(K-1)
  500          CONTINUE
               IF (WANTV) CALL DROT(P,V(1,K),1,V(1,M),1,CS,SN)
  510       CONTINUE
         GO TO 610
C
C        SPLIT AT NEGLIGIBLE S(L).
C
  520    CONTINUE
            F = E(L-1)
            E(L-1) = 0.0D0
            DO 530 K = L, M
               T1 = S(K)
               CALL DROTGX(T1,F,CS,SN)
               S(K) = T1
               F = -SN*E(K)
               E(K) = CS*E(K)
               IF (WANTU) CALL DROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
  530       CONTINUE
         GO TO 610
C
C        PERFORM ONE QR STEP.
C
  540    CONTINUE
C
C           CALCULATE THE SHIFT.
C
            SCALE = DMAX1(DABS(S(M)),DABS(S(M-1)),DABS(E(M-1)),
     *                    DABS(S(L)),DABS(E(L)))
            SM = S(M)/SCALE
            SMM1 = S(M-1)/SCALE
            EMM1 = E(M-1)/SCALE
            SL = S(L)/SCALE
            EL = E(L)/SCALE
            B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0
            C = (SM*EMM1)**2
            SHIFT = 0.0D0
            IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 550
               SHIFT = DSQRT(B**2+C)
               IF (B .LT. 0.0D0) SHIFT = -SHIFT
               SHIFT = C/(B + SHIFT)
  550       CONTINUE
            F = (SL + SM)*(SL - SM) - SHIFT
            G = SL*EL
C
C           CHASE ZEROS.
C
            MM1 = M - 1
            DO 560 K = L, MM1
               CALL DROTGX(F,G,CS,SN)
               IF (K .NE. L) E(K-1) = F
               F = CS*S(K) + SN*E(K)
               E(K) = CS*E(K) - SN*S(K)
               G = SN*S(K+1)
               S(K+1) = CS*S(K+1)
               IF (WANTV) CALL DROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
               CALL DROTGX(F,G,CS,SN)
               S(K) = F
               F = CS*E(K) + SN*S(K+1)
               S(K+1) = -SN*E(K) + CS*S(K+1)
               G = SN*E(K+1)
               E(K+1) = CS*E(K+1)
               IF (WANTU .AND. K .LT. N)
     *            CALL DROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
  560       CONTINUE
            E(M-1) = F
            ITER = ITER + 1
         GO TO 610
C
C        CONVERGENCE.
C
  570    CONTINUE
C
C           MAKE THE SINGULAR VALUE  POSITIVE.
C
            IF (S(L) .GE. 0.0D0) GO TO 580
               S(L) = -S(L)
               IF (WANTV) CALL DSCALX(P,-1.0D0,V(1,L),1)
  580       CONTINUE
C
C           ORDER THE SINGULAR VALUE.
C
  590       IF (L .EQ. MM) GO TO 600
C           ...EXIT
               IF (S(L) .GE. S(L+1)) GO TO 600
               T = S(L)
               S(L) = S(L+1)
               S(L+1) = T
               IF (WANTV .AND. L .LT. P)
     *            CALL DSWAPX(P,V(1,L),1,V(1,L+1),1)
               IF (WANTU .AND. L .LT. N)
     *            CALL DSWAPX(N,U(1,L),1,U(1,L+1),1)
               L = L + 1
            GO TO 590
  600       CONTINUE
            ITER = 0
            M = M - 1
  610    CONTINUE
      GO TO 360
  620 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB)
      INTEGER LDX,N,P,JOB
      INTEGER JPVT(1)
      DOUBLE PRECISION X(LDX,1),QRAUX(1),WORK(1)
C
C     DQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR
C     FACTORIZATION OF AN N BY P MATRIX X.  COLUMN PIVOTING
C     BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE
C     PERFORMED AT THE USERS OPTION.
C
C     ON ENTRY
C
C        X       DOUBLE PRECISION(LDX,P), WHERE LDX .GE. N.
C                X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE
C                COMPUTED.
C
C        LDX     INTEGER.
C                LDX IS THE LEADING DIMENSION OF THE ARRAY X.
C
C        N       INTEGER.
C                N IS THE NUMBER OF ROWS OF THE MATRIX X.
C
C        P       INTEGER.
C                P IS THE NUMBER OF COLUMNS OF THE MATRIX X.
C
C        JPVT    INTEGER(P).
C                JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION
C                OF THE PIVOT COLUMNS.  THE K-TH COLUMN X(K) OF X
C                IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE
C                VALUE OF JPVT(K).
C
C                   IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL
C                                      COLUMN.
C
C                   IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN.
C
C                   IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN.
C
C                BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS
C                ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL
C                COLUMNS TO THE END.  BOTH INITIAL AND FINAL COLUMNS
C                ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY
C                FREE COLUMNS ARE MOVED.  AT THE K-TH STAGE OF THE
C                REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN
C                IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST
C                REDUCED NORM.  JPVT IS NOT REFERENCED IF
C                JOB .EQ. 0.
C
C        WORK    DOUBLE PRECISION(P).
C                WORK IS A WORK ARRAY.  WORK IS NOT REFERENCED IF
C                JOB .EQ. 0.
C
C        JOB     INTEGER.
C                JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING.
C                IF JOB .EQ. 0, NO PIVOTING IS DONE.
C                IF JOB .NE. 0, PIVOTING IS DONE.
C
C     ON RETURN
C
C        X       X CONTAINS IN ITS UPPER TRIANGLE THE UPPER
C                TRIANGULAR MATRIX R OF THE QR FACTORIZATION.
C                BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM
C                WHICH THE ORTHOGONAL PART OF THE DECOMPOSITION
C                CAN BE RECOVERED.  NOTE THAT IF PIVOTING HAS
C                BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT
C                OF THE ORIGINAL MATRIX X BUT THAT OF X
C                WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT.
C
C        QRAUX   DOUBLE PRECISION(P).
C                QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER
C                THE ORTHOGONAL PART OF THE DECOMPOSITION.
C
C        JPVT    JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE
C                ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO
C                THE K-TH COLUMN, IF PIVOTING WAS REQUESTED.
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C     DQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
C
C     BLAS DAXPYX,DDOTX,DSCALX,DSWAPX,DNRM2X
C     FORTRAN DABS,DMAX1,MIN0,DSQRT
C
C     INTERNAL VARIABLES
C
      INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU
      DOUBLE PRECISION MAXNRM,DNRM2X,TT
      DOUBLE PRECISION DDOTX,NRMXL,T
      LOGICAL NEGJ,SWAPJ
C
C
      PL = 1
      PU = 0
      IF (JOB .EQ. 0) GO TO 60
C
C        PIVOTING HAS BEEN REQUESTED.  REARRANGE THE COLUMNS
C        ACCORDING TO JPVT.
C
         DO 20 J = 1, P
            SWAPJ = JPVT(J) .GT. 0
            NEGJ = JPVT(J) .LT. 0
            JPVT(J) = J
            IF (NEGJ) JPVT(J) = -J
            IF (.NOT.SWAPJ) GO TO 10
               IF (J .NE. PL) CALL DSWAPX(N,X(1,PL),1,X(1,J),1)
               JPVT(J) = JPVT(PL)
               JPVT(PL) = J
               PL = PL + 1
   10       CONTINUE
   20    CONTINUE
         PU = P
         DO 50 JJ = 1, P
            J = P - JJ + 1
            IF (JPVT(J) .GE. 0) GO TO 40
               JPVT(J) = -JPVT(J)
               IF (J .EQ. PU) GO TO 30
                  CALL DSWAPX(N,X(1,PU),1,X(1,J),1)
                  JP = JPVT(PU)
                  JPVT(PU) = JPVT(J)
                  JPVT(J) = JP
   30          CONTINUE
               PU = PU - 1
   40       CONTINUE
   50    CONTINUE
   60 CONTINUE
C
C     COMPUTE THE NORMS OF THE FREE COLUMNS.
C
      IF (PU .LT. PL) GO TO 80
      DO 70 J = PL, PU
         QRAUX(J) = DNRM2X(N,X(1,J),1)
         WORK(J) = QRAUX(J)
   70 CONTINUE
   80 CONTINUE
C
C     PERFORM THE HOUSEHOLDER REDUCTION OF X.
C
      LUP = MIN0(N,P)
      DO 200 L = 1, LUP
         IF (L .LT. PL .OR. L .GE. PU) GO TO 120
C
C           LOCATE THE COLUMN OF LARGEST NORM AND BRING IT
C           INTO THE PIVOT POSITION.
C
            MAXNRM = 0.0D0
            MAXJ = L
            DO 100 J = L, PU
               IF (QRAUX(J) .LE. MAXNRM) GO TO 90
                  MAXNRM = QRAUX(J)
                  MAXJ = J
   90          CONTINUE
  100       CONTINUE
            IF (MAXJ .EQ. L) GO TO 110
               CALL DSWAPX(N,X(1,L),1,X(1,MAXJ),1)
               QRAUX(MAXJ) = QRAUX(L)
               WORK(MAXJ) = WORK(L)
               JP = JPVT(MAXJ)
               JPVT(MAXJ) = JPVT(L)
               JPVT(L) = JP
  110       CONTINUE
  120    CONTINUE
         QRAUX(L) = 0.0D0
         IF (L .EQ. N) GO TO 190
C
C           COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L.
C
            NRMXL = DNRM2X(N-L+1,X(L,L),1)
            IF (NRMXL .EQ. 0.0D0) GO TO 180
               IF (X(L,L) .NE. 0.0D0) NRMXL = DSIGN(NRMXL,X(L,L))
               CALL DSCALX(N-L+1,1.0D0/NRMXL,X(L,L),1)
               X(L,L) = 1.0D0 + X(L,L)
C
C              APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS,
C              UPDATING THE NORMS.
C
               LP1 = L + 1
               IF (P .LT. LP1) GO TO 170
               DO 160 J = LP1, P
                  T = -DDOTX(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
                  CALL DAXPYX(N-L+1,T,X(L,L),1,X(L,J),1)
                  IF (J .LT. PL .OR. J .GT. PU) GO TO 150
                  IF (QRAUX(J) .EQ. 0.0D0) GO TO 150
                     TT = 1.0D0 - (DABS(X(L,J))/QRAUX(J))**2
                     TT = DMAX1(TT,0.0D0)
                     T = TT
                     TT = 1.0D0 + 0.05D0*TT*(QRAUX(J)/WORK(J))**2
                     IF (TT .EQ. 1.0D0) GO TO 130
                        QRAUX(J) = QRAUX(J)*DSQRT(T)
                     GO TO 140
  130                CONTINUE
                        QRAUX(J) = DNRM2X(N-L,X(L+1,J),1)
                        WORK(J) = QRAUX(J)
  140                CONTINUE
  150             CONTINUE
  160          CONTINUE
  170          CONTINUE
C
C              SAVE THE TRANSFORMATION.
C
               QRAUX(L) = X(L,L)
               X(L,L) = -NRMXL
  180       CONTINUE
  190    CONTINUE
  200 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DAXPYX(N,DA,DX,INCX,DY,INCY)
C
C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DA
      INTEGER I,INCX,INCY,IXIY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF (DA .EQ. 0.0D0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DY(IY) + DA*DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,4)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DY(I) + DA*DX(I)
   30 CONTINUE
      IF( N .LT. 4 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
        DY(I) = DY(I) + DA*DX(I)
        DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
        DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
        DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
   50 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION DDOTX(N,DX,INCX,DY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DTEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      DDOTX = 0.0D0
      DTEMP = 0.0D0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DTEMP = DTEMP + DX(IX)*DY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      DDOTX = DTEMP
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DTEMP + DX(I)*DY(I)
   30 CONTINUE
      IF( N .LT. 5 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
     *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
   50 CONTINUE
   60 DDOTX = DTEMP
      RETURN
      END
C---------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION DNRM2X ( N, DX, INCX)
      INTEGER          NEXT
      DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
      DATA   ZERO, ONE /0.0D0, 1.0D0/
C
C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
C     INCREMENT INCX .
C     IF    N .LE. 0 RETURN WITH RESULT = 0.
C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C           C.L.LAWSON, 1978 JAN 08
C
C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
C     HOPEFULLY APPLICABLE TO ALL MACHINES.
C         CUTLO = MAXIMUM OF  DSQRT(U/EPS)  OVER ALL KNOWN MACHINES.
C         CUTHI = MINIMUM OF  DSQRT(V)      OVER ALL KNOWN MACHINES.
C     WHERE
C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
C
C     BRIEF OUTLINE OF ALGORITHM..
C
C     PHASE 1    SCANS ZERO COMPONENTS.
C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C     VALUES FOR CUTLO AND CUTHI..
C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
C                   UNIVAC AND DEC AT 2**(-103)
C                   THUS CUTLO = 2**(-51) = 4.44089E-16
C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C                   THUS CUTHI = 2**(63.5) = 1.30438E19
C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
      DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C
      IF(N .GT. 0) GO TO 10
         DNRM2X  = ZERO
         GO TO 300
C
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
      NN = N * INCX
C                                                 BEGIN MAIN LOOP
      I = 1
   20    GO TO NEXT,(30, 50, 70, 110)
   30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
      ASSIGN 50 TO NEXT
      XMAX = ZERO
C
C                        PHASE 1.  SUM IS ZERO
C
   50 IF( DX(I) .EQ. ZERO) GO TO 200
      IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
C
C                                PREPARE FOR PHASE 2.
      ASSIGN 70 TO NEXT
      GO TO 105
C
C                                PREPARE FOR PHASE 4.
C
  100 I = J
      ASSIGN 110 TO NEXT
      SUM = (SUM / DX(I)) / DX(I)
  105 XMAX = DABS(DX(I))
      GO TO 115
C
C                   PHASE 2.  SUM IS SMALL.
C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
   70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
C
C                     COMMON CODE FOR PHASES 2 AND 4.
C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
C
  110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
         SUM = ONE + SUM * (XMAX / DX(I))**2
         XMAX = DABS(DX(I))
         GO TO 200
C
  115 SUM = SUM + (DX(I)/XMAX)**2
      GO TO 200
C
C
C                  PREPARE FOR PHASE 3.
C
   75 SUM = (SUM * XMAX) * XMAX
C
C
C     FOR REAL OR D.P. SET HITEST = CUTHI/N
C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
C
   85 HITEST = CUTHI/FLOAT( N )
C
C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
C
      DO 95 J =I,NN,INCX
      IF(DABS(DX(J)) .GE. HITEST) GO TO 100
   95    SUM = SUM + DX(J)**2
      DNRM2X = DSQRT( SUM )
      GO TO 300
C
  200 CONTINUE
      I = I + INCX
      IF ( I .LE. NN ) GO TO 20
C
C              END OF MAIN LOOP.
C
C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
      DNRM2X = XMAX * DSQRT(SUM)
  300 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DROTGX(DA,DB,C,S)
C
C     CONSTRUCT GIVENS PLANE ROTATION.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z
C
      ROE = DB
      IF( DABS(DA) .GT. DABS(DB) ) ROE = DA
      SCALE = DABS(DA) + DABS(DB)
      IF( SCALE .NE. 0.0D0 ) GO TO 10
         C = 1.0D0
         S = 0.0D0
         R = 0.0D0
         GO TO 20
   10 R = SCALE*DSQRT((DA/SCALE)**2 + (DB/SCALE)**2)
      R = DSIGN(1.0D0,ROE)*R
      C = DA/R
      S = DB/R
   20 Z = 1.0D0
      IF( DABS(DA) .GT. DABS(DB) ) Z = S
      IF( DABS(DB) .GE. DABS(DA) .AND. C .NE. 0.0D0 ) Z = 1.0D0/C
      DA = R
      DB = Z
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE  DSCALX(N,DA,DX,INCX)
C
C     SCALES A VECTOR BY A CONSTANT.
C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DA,DX(1)
      INTEGER I,INCX,M,MP1,N,NINCX
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I = 1,NINCX,INCX
        DX(I) = DA*DX(I)
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DX(I) = DA*DX(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DX(I) = DA*DX(I)
        DX(I + 1) = DA*DX(I + 1)
        DX(I + 2) = DA*DX(I + 2)
        DX(I + 3) = DA*DX(I + 3)
        DX(I + 4) = DA*DX(I + 4)
   50 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE  DSWAPX (N,DX,INCX,DY,INCY)
C
C     INTERCHANGES TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DTEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C         TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DTEMP = DX(IX)
        DX(IX) = DY(IY)
        DY(IY) = DTEMP
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C       CLEAN-UP LOOP
C
   20 M = MOD(N,3)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DX(I)
        DX(I) = DY(I)
        DY(I) = DTEMP
   30 CONTINUE
      IF( N .LT. 3 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,3
        DTEMP = DX(I)
        DX(I) = DY(I)
        DY(I) = DTEMP
        DTEMP = DX(I + 1)
        DX(I + 1) = DY(I + 1)
        DY(I + 1) = DTEMP
        DTEMP = DX(I + 2)
        DX(I + 2) = DY(I + 2)
        DY(I + 2) = DTEMP
   50 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
C
C     SUBROUTINES FOR UPDATING THE QR DECOMPOSITION OF A MATRIX.
C
C
C LET A=Q*R BE A QR DECOMPOSITION OF AN M BY N MATRIX A, M>=N, I.E.
C Q IS M BY N AND HAS ORTHONORMAL COLUMNS, AND R IS N BY N AND RIGHT
C TRIANGULAR. THEN
C
C   DDELC  UPDATES Q AND R WHEN A COLUMN OF MATRIX A IS DELETED.
C
C   DDELR  UPDATES Q AND R WHEN A ROW OF MATRIX A IS DELETED. FOR
C          THIS SUBROUTINE M>N IS NECESSARY.
C
C   DINSC  UPDATES Q AND R WHEN A COLUMN IS INSERTED INTO MATRIX A.
C          FOR THIS SUBROUTINE M>N IS NECESSARY.
C
C   DINSR  UPDATES Q AND R WHEN A ROW IS INSERTED INTO MATRIX A.
C
C   DRNK1  UPDATES Q AND R WHEN A RANK ONE MATRIX IS ADDED TO MATRIX A.
C
C   DRRPM  UPDATES Q AND R WHEN COLUMNS OF MATRICES A AND R ARE
C          PERMUTED. THIS SUBROUTINE CAN BE USED TO DETERMINE THE RANK
C          OF MATRIX A AND TO SOLVE THE SUBSET SELECTIONPROBLEM.
C
C THE FIRST FIVE OF THE ABOVE SUBROUTINES ARE IMPLEMENTATIONS OF
C SLIGHT MODIFICATIONS OF ALGORITHMS DESCRIBED BY J.W. DANIEL ET AL.,
C MATH. COMPUT. 30 (1976), PP. 772-795. THE SUBROUTINE DRRPM IS AN
C IMPLEMENTATION OF THE RANK REVEALING QR DECOMPOSITION SCHEME
C DESCRIBED BY T.F. CHAN, LIN. ALG. APPL., 88/89 (1987), PP. 67-82.
C
C
C AUXILIARY SUBROUTINES
C
C   DINVIT  FOR INVERSE ITERATION WITH R'*R, WHERE THE PRIME DENOTES
C           TRANSPOSITION.
C
C   DORTHO  FOR A VECTOR V THIS SUBROUTINE COMPUTES Q'*V AND
C           (I-Q*Q')*V WITH REORTHOGONALIZATION.
C
C   DORTHX  FOR AN AXIS VECTOR V THIS SUBROUTINE COMPUTES (I-Q*Q')*V
C           WITH REORTHOGONALIZATION.
C
C   DTRLSL  SOLVES LOWER TRIANGULAR SYSTEM OF EQUATIONS WITH
C           FREQUENT RESCALINGS IN ORDER TO AVOID OVERFLOW.
C
C   DTRUSL  SOLVES UPPER TRIANGULAR SYSTEM OF EQUATIONS WITH
C           FREQUENT RESCALINGS IN ORDER TO AVOID OVERFLOW.
C
C
C BLAS WRITTEN TO VECTORIZE WELL ON AN IBM 3090-200VF COMPUTER WHEN
C COMPILED WITH A VS FORTRAN 2.3.0 COMPILER. COMMENTS THROUGHOUT THE
C CODE CONCERNING VECTORIZATION REFER TO COMPILATION BY THIS COMPILER
C WITHOUT USING SPECIAL VECTORIZATION DIRECTIVES.
C
C   DADFGR, DAPLRC, DAPLRR, DAPLRV, DAPX, DASUM, DATPX, DAXPY, DCOPY,
C   DDOT, DMINRW, DSCAL, DSHFTD, DSHFTU, DZERO.
C
C
C BLAS FROM LINPACK APPENDED
C
C   DNRM2.
C
C
C THIS VERSION DATED 02/02/88. PLEASE SEND COMMENTS OR SUGGESTIONS TO
C
C                L. REICHEL
C                BERGEN SCIENTIFIC CENTRE
C                ALLEGATEN 36
C                N-5007 BERGEN, NORWAY
C
C   E-MAIL: NA.REICHEL@SCORE.STANFORD.EDU OR FSCLR@NOBERGEN.BITNET
C
C---------------------------------------------------------------------
      SUBROUTINE DDELC(Q,LDQ,M,N,R,LDR,K,V,INFO)
      INTEGER LDQ,M,N,LDR,K,INFO
      DOUBLE PRECISION Q(LDQ,N),R(LDR,N),V(M)
C
C     LET THE MATRICES Q AND R BE A QR DECOMPOSITION OF A:=Q*R.
C     GIVEN Q AND R, DDELC UPDATES THIS QR DECOMPOSITION WHEN THE
C     KTH COLUMN OF MATRIX A IS DELETED.
C
C     ON ENTRY
C
C        Q       REAL*8 (LDQ,N)
C                THE COLUMNS OF Q ARE ASSUMED ORTHONORMAL.
C
C        LDQ     INTEGER
C                LDQ IS THE LEADING DIMENSION OF THE ARRAY Q.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX Q.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF THE MATRICES Q AND R,
C                AND THE NUMBER OF ROWS OF THE MATRIX R. N<=M REQUIRED.
C
C        R       REAL*8 (LDR,N)
C                R IS AN UPPER TRIANGULAR N BY N MATRIX. A:=Q*R IS THE
C                MATRIX WHOSE KTH COLUMN IS DELETED.
C
C        LDR     INTEGER
C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
C
C        K       INTEGER
C                COLUMN K OF MATRIX A:=Q*R IS DELETED.
C
C        INFO    INTEGER
C
C                INFO=1    V CONTAINS THE DELETED COLUMN ON RETURN.
C
C                INFO=0    THE CONTENT OF V ON RETURN IS THE SAME AS
C                          ON ENTRY. V IS NOT REFERRED TO DURING THE
C                          COMPUTATIONS.
C
C     ON RETURN
C
C        Q       Q IS THE M BY N-1 MATRIX WITH ORTHONORMAL COLUMNS IN
C                A QR DECOMPOSITION OF THE MATRIX OBTAINED BY DELETING
C                COLUMN K OF MATRIX A.
C
C        R       R IS THE UPPER TRIANGULAR N-1 BY N-1 MATRIX IN A QR
C                DECOMPOSITION OF THE MATRIX OBTAINED BY DELETING
C                COLUMN K OF MATRIX A.
C
C        V       REAL*8 (N)
C                V CONTAINS THE DELETED COLUMN IF INFO=1 ON ENTRY.
C
C        INFO    INFO=0     SUCCESSFUL TERMINATION.
C
C                INFO=9     ERROR EXIT. K<1 OR K>N OR M<N OR LDQ<M
C                           OR LDR<N. NO COMPUTATIONS HAVE BEEN
C                           CARRIED OUT. ON RETURN Q AND R ARE THE
C                           SAME AS ON ENTRY.
C
C     DDELC USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     BLAS: DADFGR,DAPLRC,DAPLRR,DAXPY,DCOPY,DZERO
C
C     INTERNAL VARIABLES
C
      INTEGER L,N1
      DOUBLE PRECISION CS,SN
C
C     ERROR EXIT. IF K<1 OR K>N OR M<N OR LDQ<M OR LDR<N, SET FLAG AND
C     RETURN.
C
      IF(K.LT.1.OR.K.GT.N.OR.M.LT.N.OR.LDQ.LT.M.OR.LDR.LT.N)THEN
        INFO=9
        RETURN
      ENDIF
C
C     COMPUTE V IF INFO=1.
C
      IF(INFO.EQ.1)THEN
        CALL DZERO(V,M)
        DO 10 L=1,K
          CALL DAXPY(Q(1,L),M,R(L,K),V)
   10   CONTINUE
      ENDIF
C
C     UPDATE R AND Q.
C
      INFO=0
      N1=N-1
      DO 20 L=K,N1
        CALL DADFGR(R(L,L+1),R(L+1,L+1),CS,SN)
        CALL DAPLRR(CS,SN,R,LDR,N,L,L+2)
        CALL DAPLRC(CS,SN,M,Q(1,L),Q(1,L+1))
   20 CONTINUE
C
C     MOVE COLUMNS OF R.
C
      DO 30 L=K,N1
        CALL DCOPY(R(1,L+1),L,R(1,L))
   30 CONTINUE
C
C     SET NTH COLUMNS TO ZERO.
C
      CALL DZERO(R(1,N),N)
      CALL DZERO(Q(1,N),M)
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DDELR(Q,LDQ,M,N,R,LDR,K,U,WK,INFO)
      INTEGER K,LDQ,M,N,LDR,INFO
      DOUBLE PRECISION Q(LDQ,N),R(LDR,N),U(N),WK(2*M+N)
C
C     LET THE MATRICES Q AND R BE A QR DECOMPOSITION OF A:=Q*R.
C     GIVEN Q AND R, DDELR UPDATES THIS QR DECOMPOSITION WHEN THE
C     KTH ROW OF MATRIX A IS DELETED.
C
C     ON ENTRY
C
C        Q       REAL*8 (LDQ,N)
C                THE COLUMNS OF Q ARE ASSUMED ORTHONORMAL.
C
C        LDQ     INTEGER
C                LDQ IS THE LEADING DIMENSION OF THE ARRAY Q.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX Q.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF MATRIX Q, AND
C                THE NUMBER OF ROWS AND COLUMNS OF MATRIX R.
C                N<M REQUIRED.
C
C        R       REAL*8 (LDR,N)
C                R IS AN UPPER TRIANGULAR N BY N MATRIX. A:=Q*R IS
C                THE MATRIX WHOSE KTH ROW IS DELETED.
C
C        LDR     INTEGER
C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
C
C        K       INTEGER
C                ROW K OF MATRIX A:=Q*R IS DELETED.
C
C        WK      REAL*8 (2*M+N)
C                WK IS A SCRATCH ARRAY.
C
C     ON RETURN
C
C        Q       Q IS THE M-1 BY N MATRIX WITH ORTHONORMAL COLUMNS IN
C                A QR DECOMPOSITION OF THE MATRIX OBTAINED BY DELETING
C                ROW K OF MATRIX A.
C
C        R       R IS THE UPPER TRIANGULAR N BY N MATRIX IN A QR
C                DECOMPOSITION OF THE MATRIX OBTAINED BY DELETING
C                ROW K OF MATRIX A.
C
C        U       REAL*8 (N)
C                U(1:N) CONTAINS THE DELETED ROW.
C
C        INFO    INTEGER
C
C                INFO=0    SUCCESSFUL TERMINATION.
C
C                INFO=8    ERROR EXIT. INFO IS SET BY SUBROUTINE
C                          DORTHX, WHICH WAS UNABLE TO DETERMINE A
C                          VECTOR ORTHOGONAL TO THE COLUMNS OF Q. ON
C                          RETURN Q AND R ARE THE SAME AS ON ENTRY.
C
C                INFO=9    ERROR EXIT. K<1 OR K>M OR M<=N OR LDQ<M OR
C                          LDR<N. NO COMPUTATIONS HAVE BEEN CARRIED
C                          OUT. ON RETURN Q AND R ARE THE SAME AS ON
C                          ENTRY.
C
C     DDELR USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     BLAS: DADFGR,DAPLRC,DAPLRV,DSHFTD,DMINRW,DSCAL,DZERO
C     OTHER: DORTHX
C
C     INTERNAL VARIABLES
C
      INTEGER L,MINRW,M1
      DOUBLE PRECISION CS,SN,RHO
C
C     ERROR EXIT. IF K<1 OR K>M OR M<=N OR LDQ<M OR LDR<N, SET FLAG
C     AND RETURN.
C
      IF(K.LT.1.OR.K.GT.M.OR.M.LE.N.OR.LDQ.LT.M.OR.LDR.LT.N)THEN
        INFO=9
        RETURN
      ENDIF
C
C     ORTHOGONALIZE THE KTH AXIS VECTOR AGAINST THE COLUMNS OF Q. THE
C     SO OBTAINED VECTOR IS STORED IN WK(1:M). WK(M+1:2*M+N) IS A
C     SCRATCH ARRAY FOR SUBROUTINE DORTHX.
C
      CALL DORTHX(Q,LDQ,M,N,K,RHO,WK,WK(M+1),INFO)
C
      IF(INFO.NE.0)THEN
C
C       IF INFO>0 THEN THE KTH AXIS VECTOR WAS IN SPAN(Q) NUMERICALLY.
C       DETERMINE ROW OF Q OF LEAST LENGTH.
C
        CALL DMINRW(Q,LDQ,M,N,WK,MINRW)
        CALL DORTHX(Q,LDQ,M,N,MINRW,RHO,WK,WK(M+1),INFO)
C
C       ERROR EXIT. IF INFO>0 THEN NO VECTOR WK(1:M) SUCH THAT Q'*WK=0
C       HAS BEEN FOUND.
C
        IF(INFO.NE.0)RETURN
        RHO=0D0
      ENDIF
C
C     LOOP IS ELIGIBLE FOR VECTORIZATION, BUT DOES NOT VECTORIZE
C     WITHOUT SPECIAL COMPILER DIRECTIVES.
C
      DO 10 L=1,N
        U(L)=Q(K,L)
   10 CONTINUE
      CALL DSHFTD(WK(K),M-K+1)
      M1=M-1
      DO 20 L=N,1,-1
        CALL DSHFTD(Q(K,L),M-K+1)
C
C       INTRODUCE ZERO COMPONENTS IN ROW OF Q BY APPLYING GIVENS
C       REFLECTORS.
C
        CALL DADFGR(RHO,U(L),CS,SN)
        CALL DAPLRV(CS,SN,U(L),1,R(L,L),LDR,N-L+1)
        CALL DAPLRC(CS,SN,M1,WK,Q(1,L))
        Q(M,L)=0D0
   20 CONTINUE
      CALL DSCAL(U,N,RHO)
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DINSC(Q,LDQ,M,N,R,LDR,K,V,RCOND,WK,INFO)
      INTEGER K,LDQ,M,N,LDR,INFO
      DOUBLE PRECISION Q(LDQ,N),R(LDR,N),V(M),WK(2*M+2*N),RCOND
C
C     LET THE MATRICES Q AND R BE A QR DECOMPOSITION OF A:=Q*R.
C     GIVEN Q AND R, DINSC UPDATES THIS QR DECOMPOSITION WHEN THE
C     VECTOR V IS INSERTED BETWEEN COLUMNS K-1 AND K OF MATRIX A.
C
C     ON ENTRY
C
C        Q       REAL*8 (LDQ,N)
C                THE COLUMNS OF Q ARE ASSUMED ORTHONORMAL.
C
C        LDQ     INTEGER
C                LDQ IS THE LEADING DIMENSION OF THE ARRAY Q.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX Q.
C
C        N       INTEGER
C                N-1 IS THE NUMBER OF COLUMNS OF THE MATRICES Q AND R,
C                AND THE NUMBER OF ROWS AND COLUMNS OF MATRIX R.
C                N<=M REQUIRED.
C
C        R       REAL*8 (LDR,N)
C                R IS AN UPPER TRIANGULAR N-1 BY N-1 MATRIX. A:=Q*R IS
C                THE MATRIX TO BE UPDATED.
C
C        LDR     INTEGER
C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
C
C        K       INTEGER
C                VECTOR V IS INSERTED BETWEEN COLUMNS K-1 AND K OF
C                MATRIX A:=Q*R, 1<=K<=N.
C
C        V       REAL*8 (M)
C                VECTOR V IS INSERTED BETWEEN COLUMNS K-1 AND K OF
C                MATRIX A:=Q*R.
C
C        RCOND   REAL*8
C                DO NOT INSERT V IF THE RECIPROCAL CONDITION NUMBER OF
C                THE MATRIX OBTAINED BY APPENDING COLUMN V/LENGHT(V) TO
C                MATRIX Q IS LESS THAN RCOND.
C
C        WK      REAL*8 (2*M+2*N)
C                WK IS A SCRATCH ARRAY.
C
C     ON RETURN
C
C        Q       Q IS THE M BY N MATRIX WITH ORTHONORMAL COLUMNS IN A
C                QR DECOMPOSITION OF THE MATRIX OBTAINED BY INSERTING
C                COLUMN V INTO MATRIX A.
C
C        R       R IS THE UPPER TRIANGULAR N BY N MATRIX IN A QR
C                DECOMPOSITION OF THE MATRIX OBTAINED BY INSERTING
C                COLUMN V INTO MATRIX A.
C
C        RCOND   RCOND IS THE RECIPROCAL CONDITION NUMBER OF THE
C                M BY N MATRIX WHOSE COLUMNS ARE V/LENGHT(V) AND THOSE
C                OF Q WHERE V AND Q ARE AS DEFINED ON ENTRY.
C
C        INFO    INTEGER
C
C                INFO=0    SUCCESSFUL TERMINATION.
C
C                INFO=1    ERROR EXIT. V WAS FOUND TO LIE IN THE RANGE
C                          OF Q NUMERICALLY. SUBROUTINE DORTHO WAS
C                          UNABLE TO DETERMINE A VECTOR ORTHOGONAL TO
C                          THE COLUMNS OF Q. Q, R AND V ARE THE SAME
C                          AS ON ENTRY.
C
C                INFO=2    THE M BY N MATRIX WITH COLUMNS Q AND
C                          V/LENGTH(V) WITH Q AND V AS GIVEN ON ENTRY
C                          HAS RECIPROCAL CONDITION NUMBER SMALLER
C                          THAN RCOND AS GIVEN ON ENTRY. Q, R AND V
C                          ARE THE SAME AS ON ENTRY.
C
C                INFO=9    ERROR EXIT. K>0 OR K<N OR M<N OR LDQ<M OR
C                          LDR<N. Q, R, AND V ARE THE SAME AS ON ENTRY.
C
C     DINSC USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     BLAS: DADFGR,DAPLRC,DAPLRR,DCOPY
C     OTHER: DORTHO
C
C     INTERNAL VARIABLES
C
      INTEGER L,N1
      DOUBLE PRECISION CS,SN
C
C     ERROR EXIT. IF K<1 OR K>N OR M<N OR LDQ<M OR LDR<N, SET FLAG AND
C     RETURN.
C
      IF(K.LT.1.OR.K.GT.N.OR.M.LT.N.OR.LDQ.LT.M.OR.LDR.LT.N)THEN
        INFO=9
        RETURN
      ENDIF
C
C     ORTHOGONALIZE V AND INSERT ORTHOGONALIZED VECTOR INTO Q.
C
      N1=N-1
      CALL DORTHO(Q,LDQ,M,N1,V,WK,WK(N+1),RCOND,WK(M+N+1),INFO)
C
C     GENERALLY IT IS NOT MEANINGFUL TO INSERT A COLUMN THAT IS
C     A LINEAR COMBINATION OF ALREADY EXISTING COLUMNS. EXIT IF
C     INFO IS NONVANISHING.
C
      IF(INFO.NE.0)RETURN
C
C     LET VECTOR WK(N+1:M+N+1), WHICH IS ORTHOGONAL TO Q, BE A NEW
C     LAST COLUMN OF Q. MOVE COLUMNS N-1, N-2, ..., K OF R.
C
      CALL DCOPY(WK(N+1),M,Q(1,N))
      DO 10 L=N1,K,-1
        CALL DCOPY(R(1,L),L,R(1,L+1))
   10 CONTINUE
C
C     LOOP IS ELIGIBLE FOR VECTORIZATION, BUT DOES NOT VECTORIZE
C     WITHOUT SPECIAL COMPILER DIRECTIVES.
C
      DO 20 L=K+1,N
        R(L,L)=0D0
   20 CONTINUE
C
C     ELEMENTS WK(K+1:N) ARE ANNIHILATED BY GIVENS REFLECTORS SO THAT
C     R BECOMES UPPER TRIANGULAR.
C
      DO 30 L=N1,K,-1
        CALL DADFGR(WK(L),WK(L+1),CS,SN)
        CALL DAPLRR(CS,SN,R,LDR,N,L,L+1)
        CALL DAPLRC(CS,SN,M,Q(1,L),Q(1,L+1))
   30 CONTINUE
      CALL DCOPY(WK,K,R(1,K))
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DINSR(Q,LDQ,M,N,R,LDR,K,U,WK,INFO)
      INTEGER LDQ,M,N,LDR,K,INFO
      DOUBLE PRECISION Q(LDQ,N),R(LDR,N),U(N),WK(M)
C
C     LET THE MATRICES Q AND R BE A QR DECOMPOSITION OF A:=Q*R.
C     GIVEN Q AND R, DINSR UPDATES THIS QR DECOMPOSITION WHEN THE
C     VECTOR U IS INSERTED BETWEEN ROWS K-1 AND K OF MATRIX A.
C
C     ON ENTRY
C
C        Q       REAL*8 (LDQ,N)
C                THE COLUMNS OF Q ARE ASSUMED ORTHONORMAL.
C
C        LDQ     INTEGER
C                LDQ IS THE LEADING DIMENSION OF THE ARRAY Q.
C
C        M       INTEGER
C                M-1 IS THE NUMBER OF ROWS OF THE MATRIX Q.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF THE MATRICES Q AND R,
C                AND THE NUMBER OF ROWS OF R. N<M REQUIRED.
C
C        R       REAL*8 (LDR,N)
C                R IS AN UPPER TRIANGULAR N BY N MATRIX. A:=Q*R IS THE
C                MATRIX WHICH IS UPDATED.
C
C        LDR     INTEGER
C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
C
C        K       INTEGER
C                VECTOR U IS INSERTED BETWEEN ROWS K-1 AND K OF
C                MATRIX A:=Q*R, 1<=K<=M.
C
C        U       REAL*8 (N)
C                VECTOR CONTAINING ELEMENTS OF ROW TO BE INSERTED.
C
C        WK      REAL*8 (M)
C                WK IS A SCRATCH ARRAY.
C
C     ON RETURN
C
C        Q       Q IS THE M BY N MATRIX WITH ORTHONORMAL COLUMNS IN A
C                QR DECOMPOSITION OF THE MATRIX OBTAINED BY INSERTING
C                ROW U INTO MATRIX A.
C
C        R       R IS THE UPPER TRIANGULAR N BY N MATRIX IN A QR
C                DECOMPOSITION OF THE MATRIX OBTAINED BY INSERTING
C                ROW U INTO MATRIX A.
C
C        INFO    INTEGER
C
C                INFO=0     SUCCESSFUL TERMINATION.
C
C                INFO=9     ERROR EXIT. M<=N OR K<1 OR K>M OR LDQ<M
C                           OR LDR<N.  NO COMPUTATIONS HAVE BEEN
C                           CARRIED OUT. ON RETURN Q, R AND U ARE THE
C                           SAME AS ON ENTRY.
C
C     DINSR USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     BLAS: DADFGR,DAPLRC,DAPLRV,DSHFTU,DZERO
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION CS,SN
      INTEGER L,N1
C
C     ERROR EXIT. IF K<1 OR K>M OR M<=N OR LDQ<M OR LDR<N, SET FLAG
C     AND RETURN.
C
      IF(K.LT.1.OR.K.GT.M.OR.M.LE.N.OR.LDQ.LT.M.OR.LDR.LT.N)THEN
        INFO=9
        RETURN
      ENDIF
C
      INFO=0
      CALL DZERO(WK,M)
      WK(K)=1D0
C
      DO 10 L=1,N
C
C       SHIFT ELEMENTS OF VECTOR (Q(K,L),Q(K+1,L),...,Q(M-1,L))'.
C
        CALL DSHFTU(Q(K,L),M-K)
        Q(K,L)=0D0
C
C       UPDATE R AND Q.
C
        CALL DADFGR(R(L,L),U(L),CS,SN)
        CALL DAPLRV(CS,SN,R(L,L+1),LDR,U(L+1),1,N-L)
        CALL DAPLRC(CS,SN,M,Q(1,L),WK)
   10 CONTINUE
C
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DRNK1(Q,LDQ,M,N,R,LDR,U,V,WK,INFO)
      INTEGER LDQ,M,N,LDR,INFO
      DOUBLE PRECISION Q(LDQ,N),R(LDR,N),U(N),V(M),WK(2*M+2*N+1)
C
C     LET THE MATRICES Q AND R BE A QR DECOMPOSITION OF A:=Q*R.
C     GIVEN Q AND R, DRNK1 COMPUTES A QR DECOMPOSITION OF THE RANK ONE
C     MODIFICATION B:=A+V*U' OF A. HERE U AND V ARE VECTORS AND THE
C     PRIME DENOTES TRANSPOSITION.
C
C     ON ENTRY
C
C        Q       REAL*8 (LDQ,N)
C                THE COLUMNS OF Q ARE ASSUMED ORTHONORMAL.
C
C        LDQ     INTEGER
C                LDQ IS THE LEADING DIMENSION OF THE ARRAY Q.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX Q.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF THE MATRICES Q AND R,
C                AND THE NUMBER OF ROWS OF THE MATRIX R.
C                N<=M REQUIRED.
C
C        R       REAL*8 (LDR,N)
C                R IS AN UPPER TRIANGULAR N BY N MATRIX. A:=Q*R IS THE
C                MATRIX WHICH IS UPDATED BY A RANK ONE MODIFICATION.
C
C        LDR     INTEGER
C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
C
C        U       REAL*8 (N)
C                U DEFINES THE RANK ONE MODIFICATION OF A:=Q*R.
C
C        V       REAL*8 (M)
C                V DEFINES THE RANK ONE MODIFICATION OF A:=Q*R.
C
C        WK      REAL*8 (2*M+2*N+1)
C                WK IS A SCRATCH ARRAY.
C
C     ON RETURN
C
C        Q       Q IS THE MATRIX WITH ORTHONORMAL COLUMNS IN A QR
C                DECOMPOSITION OF THE MATRIX B:=A+V*U'.
C
C        R       R IS THE UPPER TRIANGULAR MATRIX IN A QR
C                DECOMPOSITION OF THE MATRIX B:=A+V*U'.
C
C        INFO    INTEGER
C
C                INFO=0    SUCCESSFUL TERMINATION. VECTOR V WAS BY
C                          SUBROUTINE DORTHO FOUND TO BE LINEARLY
C                          INDEPENDENT OF THE COLUMNS OF Q AND MADE
C                          ORTHOGONAL TO Q.
C
C                INFO=1    SUCCESSFUL TERMINATION. VECTOR V WAS BY
C                          SUBROUTINE DORTHO FOUND TO LIE IN SPAN(Q).
C                          THIS SIMPLIFIES THE UPDATING SLIGHTLY.
C
C                INFO=9    ERROR EXIT. M<N OR LDQ<M OR LDR<N. NO
C                          COMPUTATIONS HAVE BEEN CARRIED OUT. Q, R
C                          U AND V ARE THE SAME AS ON ENTRY.
C
C     DRNK1 USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     BLAS: DADFGR,DAPLRC,DAPLRR
C     OTHER: DORTHO
C
C     INTERNAL VARIABLES
C
      INTEGER K,N1
      DOUBLE PRECISION RCOND,CS,SN,TMP
C
C     ERROR EXIT. IF M<N OR LDQ<M OR LDR<N, SET FLAG AND RETURN.
C
      IF(M.LT.N.OR.LDQ.LT.M.OR.LDR.LT.N)THEN
        INFO=9
        RETURN
      ENDIF
C
      N1=N-1
      RCOND=0D0
C
C     ORTHOGONALIZE V. ON RETURN FROM SUBROUTINE DORTHO THE FOURIER
C     COEFFICIENTS OF V WITH RESPECT TO THE COLUMNS OF Q ARE IN
C     WK(1:N), THE DISTANCE BETWEEN V AND SPAN(Q) IS IN WK(N+1), AND
C     THE ORTHOGONALIZED VECTOR OBTAINED FROM V IS STORED IN
C     WK(N+2:M+N+1). WK(M+N+2:2*M+2*N+1) IS A SCRATCH ARRAY FOR
C     SUBROUTINE DORTHO.
C
      CALL DORTHO(Q,LDQ,M,N,V,WK,WK(N+2),RCOND,WK(M+N+2),INFO)
C
C     PROJECT WK(1:N+1) ONTO AXIS VECTOR E1 WITH GIVENS REFLECTORS.
C     IF INFO=1 THEN WK(N+1)=0. INFO=0 OR INFO=1.
C
      IF(INFO.EQ.0)THEN
C
C       ZERO WK(N+1) AND APPLY REFLECTOR TO R AND Q.
C
        CALL DADFGR(WK(N),WK(N+1),CS,SN)
        TMP=R(N,N)
        R(N,N)=TMP*CS
        WK(N+1)=TMP*SN
        CALL DAPLRC(CS,SN,M,Q(1,N),WK(N+2))
      ENDIF
C
C     ZERO COMPONENTS WK(N),WK(N-1)...WK(2), AND APPLY REFLECTORS
C     TO R AND Q.
C
      DO 10 K=N1,1,-1
        CALL DADFGR(WK(K),WK(K+1),CS,SN)
        CALL DAPLRR(CS,SN,R,LDR,N,K,K)
        CALL DAPLRC(CS,SN,M,Q(1,K),Q(1,K+1))
   10 CONTINUE
C
C     LOOP IS ELIGIBLE FOR VECTORIZATION, BUT DOES NOT VECTORIZE
C     WITHOUT SPECIAL COMPILER DIRECTIVES.
C
      TMP=WK(1)
      DO 20 K=1,N
        R(1,K)=R(1,K)+TMP*U(K)
   20 CONTINUE
C
C     R IS OF UPPER HESSENBERG FORM. TRIANGULARIZE R.
C
      DO 30 K=1,N1
        CALL DADFGR(R(K,K),R(K+1,K),CS,SN)
        CALL DAPLRR(CS,SN,R,LDR,N,K,K+1)
        CALL DAPLRC(CS,SN,M,Q(1,K),Q(1,K+1))
   30 CONTINUE
      IF(INFO.EQ.0)THEN
        CALL DADFGR(R(N,N),WK(N+1),CS,SN)
        CALL DAPLRC(CS,SN,M,Q(1,N),WK(N+2))
      ENDIF
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DRRPM(Q,LDQ,M,N,R,LDR,K,IP,DELTA,NMBIT,IPOS,WK,INFO)
      INTEGER LDQ,M,N,LDR,K,IP(N),NMBIT,IPOS(0:NMBIT),INFO
      DOUBLE PRECISION Q(LDQ,N),R(LDR,N),DELTA,WK(N)
C
C     LET THE MATRICES Q AND R BE A QR DECOMPOSITION OF A*P:=Q*R,
C     WHERE P IS A PERMUTATION MATRIX. DRRPM COMPUTES A RANK REVEALING
C     QR DECOMPOSITION OF A*P ACCORDING TO A METHOD DESCRIBED BY
C     T.F. CHAN, LIN. ALG. APPL., 88/89 (1987), PP. 67-82. ON RETURN
C     DRRPM YIELDS A QR DECOMPOSITION OF A*P FOR A (DIFFERENT)
C     PERMUTATION P, SUCH THAT IF MATRIX A HAS SMALL SINGULAR VALUES,
C     THEN GENERALLY THE LOWER RIGHT HAND SIDE CORNER OF R CONTAINS
C     SMALL ELEMENTS.
C
C     ON ENTRY
C
C        Q       REAL*8 (LDQ,N)
C                THE COLUMNS OF Q ARE ASSUMED ORTHONORMAL.
C
C        LDQ     INTEGER
C                LDQ IS THE LEADING DIMENSION OF THE ARRAY Q.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX Q.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF MATRIX Q, AND
C                THE NUMBER OF ROWS AND COLUMNS OF MATRIX R.
C
C        R       REAL*8 (LDR,N)
C                R IS AN UPPER TRIANGULAR MATRIX. THE MATRICES Q
C                AND R ARE THOUGHT OF AS BEING A QR DECOMPOSITION
C                OF A MATRIX A*P, WHERE P IS A PERMUTATION MATRIX.
C
C        LDR     INTEGER
C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
C
C        K       INTEGER
C                THE ELEMENTS OF R WITH ROW AND COLUMN INDICES LARGER
C                THAN K VANISH. DRRPM WILL TYPICALLY BE CALLED FOR
C                K=N,N-1,N-2,...  UNTIL THE RANK OF MATRIX A HAS BEEN
C                FOUND.
C
C        IP      INTEGER (N)
C                IP REPRESENTS THE PERMUTATION MATRIX P. THE JTH
C                COLUMN OF A*P IS THE IP(J)TH COLUMN OF MATRIX A,
C                FOR J=1,2,...,N. IF K=N THEN GENERALLY IP(J)=J FOR
C                ALL J.
C
C        DELTA   REAL*8
C                THE COLUMNS OF R ARE NOT PERMUTED IF THE COMPUTED
C                UPPER BOUND FOR THE LEAST SINGULAR VALUE OF THE
C                LEADING PRINCIPAL K BY K SUBMATRIX OF R IS LARGER
C                THAN DELTA.
C
C        NMBIT   INTEGER
C                NUMBER OF INVERSE ITERATIONS. NMBIT HAS TO BE AT
C                LEAST ONE.
C
C        WK      REAL*8 (N)
C                WK IS A SCRATCH ARRAY.
C
C     ON RETURN
C
C        Q       Q CONTAINS ORTHONORMAL COLUMNS. MATRICES Q AND R
C                CONSTITUTE A QR DECOMPOSITION OF THE MATRIX A*P,
C                WHERE THE PERMUTATION P IS REPRESENTED BY IP. IF A
C                HAS SMALL SINGULAR VALUES, THEN R GENERALLY HAS SMALL
C                ELEMENTS IN ITS LOWER RIGHT HAND SIDE CORNER.
C
C        R       R IS THE UPPER TRIANGULAR N BY N MATRIX IN A QR
C                DECOMPOSITION OF A*P.
C
C        IP      VECTOR CONTAING PERMUTATION OF COLUMNS OF A, I.E.
C                THE QR DECOMPOSITION COMPUTED IS OF A MATRIX WHOSE
C                JTH COLUMN IS COLUMN IP(J) OF MATRIX A, J=1,2,...,N.
C
C        DELTA   DELTA IS AN UPPER BOUND FOR THE LEAST SINGULAR VALUE
C                THE K BY K LEADING PRINCIPAL SUBMATRIX OF R. DELTA IS
C                COMPUTED BY SUBROUTINE DINVIT.
C
C        IPOS    INTEGER (0:NMBIT)
C                IPOS(J) CONTAINS THE POSITION OF AN ELEMENT OF
C                LARGEST MAGNITUDE OF AN APPROXIMATE RIGHT SINGULAR
C                VECTOR CORRESPONDING TO A LEAST SINGULAR VALUE
C                OF THE LEADING K BY K SUBMATRIX OF R. IPOS(0) IS FOR
C                AN INITIAL VECTOR V OBTAINED BY A DTRCO-LIKE SCHEME.
C                IPOS(J), 1<=J<=NMBIT, IS FOR VECTORS OBTAINED BY J
C                INVERSE ITERATIONS, USING V AS A STARTING VECTOR.
C                IPOS IS SET BY SUBROUTINE DINVIT.
C
C        INFO    INTEGER
C
C        INFO    INFO=0    SUCCESSFUL TERMINATION. R IS NONSINGULAR
C                          AND NMBIT INVERSE ITERATIONS HAVE BEEN
C                          CARRIED OUT.
C
C                INFO=1    SUCCESSFUL TERMINATION. R HAS A VANISHING
C                          DIAGONAL ELEMENT. IPOS(0) CONTAINS THE
C                          POSITION OF AN ELEMENT OF LARGEST MAGNITUDE
C                          OF A RIGHT SINGULAR VECTOR CORRESPONDING TO
C                          SINGULAR VALUE 0.
C
C                INFO=2    SUCCESSFUL TERMINATION. THE COMPUTED UPPER
C                          BOUND FOR THE LEAST SINGULAR VALUE OF THE
C                          LEADING K BY K PRINCIPAL SUBMATRIX OF R
C                          IS LARGER THAN THE VALUE OF DELTA ON ENTRY.
C                          Q, R AND IP ARE THE SAME AS ON ENTRY.
C
C                INFO=8    ERROR EXIT. DUE TO UNDERFLOW ALL THE
C                          COMPONENTS OF AN APPROXIMATE RIGHT SINGULAR
C                          VECTOR VANISHED, AND COMPUTATIONS HAD TO
C                          BE DISCONTINUED. RESCALING R MAY PREVENT
C                          UNDERFLOW. THE CONTENT OF IPOS AND DELTA
C                          IS UNDEFINED.
C
C                INFO=9    ERROR EXIT. K<1 OR K>N OR M<N OR LDQ<M OR
C                          LDR<N OR NMBIT<1. NO COMPUTATIONS HAVE BEEN
C                          CARRIED OUT. Q AND R ARE THE SAME AS ON
C                          ENTRY.
C
C     DRRPM USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     BLAS: DADFGR,DAPLRC,DAPLRR,DCOPY,DZERO
C     LINPACK BLAS: DNRM2
C     OTHER: DINVIT
C
C     INTERNAL VARIABLES
C
      INTEGER I,J,K1,MAXIDX
      DOUBLE PRECISION CS,SN,TMP
C
C     ERROR EXIT. IF K<1 OR K>N OR M<N OR LDQ<M OR LDR<N OR NMBIT<1,
C     SET FLAG AND RETURN.
C
      IF(K.LT.1.OR.K.GT.N.OR.M.LT.N.OR.LDQ.LT.M.OR.LDR.LT.N.OR.
     %  NMBIT.LT.1)THEN
        INFO=9
        RETURN
      ENDIF
C
C     COMPUTE THE SINGULAR VECTOR WK CORRESPONDING TO THE LEAST
C     SINGULAR VALUE OF THE LEADING K BY K MATRIX OF R.
C
      TMP=DELTA
      CALL DINVIT(R,LDR,K,NMBIT,WK,DELTA,IPOS,INFO)
C
C     ERROR EXIT. IF INFO=8 THEN WK UNDERFLOWED IN DINVIT.
C
      IF(INFO.EQ.8)RETURN
C
C     EXIT. IF TMP<DELTA SET FLAG AND RETURN.
C
      IF(TMP.LT.DELTA)THEN
        INFO=2
        RETURN
      ENDIF
C
      IF(INFO.EQ.0)THEN
        MAXIDX=IPOS(NMBIT)
      ELSE
        MAXIDX=IPOS(0)
      ENDIF
C
C     NO REORDERING NECESSARY IF MAXIDX=K.
C
      IF(MAXIDX.EQ.K)RETURN
      K1=K-1
C
C     UPDATE PERMUTATION VECTOR IP. CYCLIC LEFT SHIFT OF
C     IP(MAXIDX),IP(MAXIDX+1),...,IP(K).
C
C     LOOP VECTORIZES.
C
      J=IP(MAXIDX)
      DO 10 I=MAXIDX,K1
        IP(I)=IP(I+1)
   10 CONTINUE
      IP(K)=J
C
C     UPDATE R. CYCLIC LEFT SHIFT OF COLUMNS MAXIDX, MAXIDX+1,...,K
C     OF R.
C
      CALL DCOPY(R(1,MAXIDX),MAXIDX,WK)
      DO 20 I=MAXIDX,K1
        CALL DCOPY(R(1,I+1),I+1,R(1,I))
   20 CONTINUE
      CALL DCOPY(WK,MAXIDX,R(1,K))
      CALL DZERO(R(MAXIDX+1,K),K-MAXIDX)
C
C     TRIANGULARIZE LEADING K BY K SUBMATRIX OF R BY APPLYING
C     GIVENS REFLECTORS. APPLY REFLECTORS TO Q.
C
      DO 30 I=MAXIDX,K1
        CALL DADFGR(R(I,I),R(I+1,I),CS,SN)
        CALL DAPLRR(CS,SN,R,LDR,N,I,I+1)
        CALL DAPLRC(CS,SN,M,Q(1,I),Q(1,I+1))
   30 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
C SUBPROGRAMS USED BY THE ABOVE SUBROUTINES.
C---------------------------------------------------------------------
      SUBROUTINE DINVIT(R,LDR,N,NMBIT,Z,DELTA,IPOS,INFO)
      INTEGER LDR,N,NMBIT,IPOS(0:NMBIT),INFO
      DOUBLE PRECISION R(LDR,N),DELTA,Z(N)
C
C     DINVIT COMPUTES AN APPROXIMATION OF THE SMALLEST SINGULAR VALUE
C     AND THE CORRESPONDING SINGULAR VECTOR OF THE RIGHT TRIANGULAR
C     MATRIX R. BY COMPUTATIONS SIMILAR TO THOSE OF SUBROUTINE DTRCO
C     OF LINPACK A VECTOR IS OBTAINED THAT GENERALLY IS A GOOD
C     APPROXIMATION OF A RIGHT SINGULAR VECTOR CORRESPONDING TO A
C     SINGULAR VALUE OF LEAST MAGNITUDE. THIS VECTOR IS IMPROVED BY
C     INVERSE ITERATIONS. FREQUENT RESCALINGS ARE CARRIED OUT IN
C     ORDER TO AVOID OVERFLOW.
C
C     ON ENTRY
C
C        R       REAL*8 (LDR,N)
C                R IS AN UPPER TRIANGULAR N BY N MATRIX.
C
C        LDR     INTEGER
C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
C
C        N       INTEGER
C                N IS THE NUMBER OF ROWS AND COLUMNS OF MATRIX R.
C
C        NMBIT   INTEGER
C                NMBIT IS THE NUMBER OF INVERSE ITERATIONS. NMBIT MUST
C                BE AT LEAST ONE. NMBIT=2 IS GENERALLY SUFFICIENT.
C
C     ON RETURN
C
C        Z       REAL*8 (N)
C                Z IS AN APPROXIMATION OF THE RIGHT SINGULAR VECTOR
C                BELONGING TO THE LEAST SINGULAR VALUE. Z IS UNSCALED.
C
C        DELTA   REAL*8
C                DELTA IS AN UPPER BOUND FOR THE LEAST SINGULAR VALUE.
C                IF Z IS A GOOD APPROXIMATION OF THE RIGHT SINGULAR
C                VECTOR BELONGING TO THE LEAST SINGULAR VALUE, THEN
C                DELTA IS A GOOD APPROXIMATION OF THE LEAST SINGULAR
C                VALUE.
C
C        IPOS    INTEGER(0:NMBIT)
C
C                IPOS(0) IS THE POSITION OF AN ELEMENT OF Z OF LARGEST
C                MAGNITUDE, WHERE Z HAS BEEN COMPUTED BY A DTRCO-LIKE
C                SCHEME (FOR INFO=0), OR WHERE Z HAS BEEN DETERMINED
C                BY OBSERVING THAT R IS SINGULAR (FOR INFO=1).
C
C                IPOS(J) IS THE POSITION OF AN ELEMENT OF Z OF LARGEST
C                MAGNITUDE WHERE Z HAS BEEN OBTAINED BY THE JTH INVERSE
C                ITERATION, J=1,2,...,NMBIT. THE ELEMENTS IPOS(J), J>0,
C                ARE UNDEFINED IF INFO IS NONVANISHING.
C
C        INFO    INTEGER
C
C        INFO    INFO=0    SUCCESSFUL TERMINATION. R IS NONSINGULAR AND
C                          NMBIT INVERSE ITERATIONS HAVE BEEN CARRIED
C                          OUT.
C
C                INFO=1    SUCCESSFUL TERMINATION. R HAS A VANISHING
C                          DIAGONAL ELEMENT. IPOS(0) CONTAINS THE
C                          POSITION OF AN ELEMENT OF LARGEST MAGNITUDE
C                          OF A RIGHT SINGULAR VECTOR CORRESPONDING TO
C                          SINGULAR VALUE 0.
C
C                INFO=8    ERROR EXIT. DUE TO UNDERFLOW ALL THE
C                          COMPONENTS OF THE APPROXIMATE RIGHT SINGULAR
C                          VECTOR Z VANISHED, AND COMPUTATIONS HAD TO
C                          BE DISCONTINUED. RESCALING R MAY PREVENT
C                          UNDERFLOW. THE CONTENT OF IPOS AND DELTA
C                          IS UNDEFINED.
C
C                INFO=9    ERROR EXIT. NMBIT WAS LESS THAN ONE ON
C                          ENTRY. NO COMPUTATIONS HAVE BEEN CARRIED
C                          OUT.
C
C     DINVIT USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     BLAS: DASUM,DCOPY,DSCAL,DTRLSL,DTRUSL,DZERO,IDMAX
C     LINPACK BLAS: DNRM2
C
C     INTERNAL VARIABLES
C
      INTEGER J,K,IDMAX
      DOUBLE PRECISION EK,S,SM,W,WK,WKM,TMP,SCALE,DASUM,DNRM2
C
C     ERROR EXIT.
C
      IF(NMBIT.LT.1)THEN
        INFO=9
        RETURN
      ENDIF
C
C     CHECK IF R(1,1) VANISHES.
C
      IF(R(1,1).EQ.0D0)THEN
        CALL DZERO(Z(2),N-1)
        Z(1)=1D0
        DELTA=0D0
        IPOS(0)=1
        INFO=1
        RETURN
      ENDIF
C
C     CHECK IF R(2,2),...,R(N,N) VANISH.
C
      DO 10 J=2,N
        IF(R(J,J).EQ.0D0)THEN
          CALL DCOPY(R(1,J),J-1,Z)
          CALL DTRUSL(R,LDR,J-1,Z,SCALE)
          Z(J)=-SCALE
          CALL DZERO(Z(J+1),N-J)
          DELTA=0D0
          IPOS(0)=IDMAX(Z,J)
          INFO=1
          RETURN
        ENDIF
   10 CONTINUE
C
C     MATRIX R IS NONSINGULAR. SOLVE TRANSPOSE(R)*Y=B WHERE THE
C     ELEMENTS OF B ARE CHOSEN TO CAUSE MAXIMAL GROWTH IN THE ELEMENTS
C     OF Y. THE CHOICE IS THE SAME AS IN SUBROUTINE DTRCO OF LINPACK.
C
      INFO=0
      EK=1D0
      CALL DZERO(Z,N)
C
      DO 20 K=1,N
        IF(Z(K).NE.0D0) EK=DSIGN(EK,-Z(K))
        IF(DABS(EK-Z(K)).GT.DABS(R(K,K)))THEN
          S=DABS(R(K,K)/(EK-Z(K)))
          CALL DSCAL(Z,N,S)
          EK=S*EK
        ENDIF
C
        WK=EK-Z(K)
        WKM=-EK-Z(K)
        S=DABS(WK)
        SM=DABS(WKM)
        WK=WK/R(K,K)
        WKM=WKM/R(K,K)
C
        IF(K.LT.N)THEN
C
C         LOOP VECTORIZES.
C
          DO 30 J=K+1,N
            SM=SM+DABS(Z(J)+WKM*R(K,J))
            Z(J)=Z(J)+WK*R(K,J)
            S=S+DABS(Z(J))
   30     CONTINUE
C
          IF(S.LT.SM)THEN
            W=WKM-WK
            WK=WKM
C
C           LOOP VECTORIZES.
C
            DO 40 J=K+1,N
              Z(J)=Z(J)+W*R(K,J)
   40       CONTINUE
          ENDIF
        ENDIF
C
        Z(K)=WK
   20 CONTINUE
C
C     NORMALIZE. DASUM IS SIMPLER THAN DNRM2 AND SUFFICES WHEN LEAST
C     SQUARES NORM IS NOT ESSENTIAL.
C
      TMP=DASUM(Z,N)
C
C     ERROR EXIT. IF LENGTH(Z)=0, SET FLAG AND RETURN.
C
      IF(TMP.EQ.0D0)THEN
        INFO=8
        RETURN
      ENDIF
C
      CALL DSCAL(Z,N,1D0/TMP)
C
C     SOLVE R*Z=Y.
C
      CALL DTRUSL(R,LDR,N,Z,SCALE)
      IPOS(0)=IDMAX(Z,N)
C
C     INVERSE ITERATION. NORMALIZE Z TO HAVE LEAST SQUARES NORM ONE IN
C     THE BEGINNING OF EACH ITERATION.
C
      DO 50 J=1,NMBIT
        TMP=DNRM2(N,Z,1)
C
C       ERROR EXIT. IF LENGTH(Z)=0, SET FLAG AND RETURN.
C
        IF(TMP.EQ.0D0)THEN
          INFO=8
          RETURN
        ENDIF
C
        CALL DSCAL(Z,N,1D0/TMP)
        CALL DTRLSL(R,LDR,N,Z,SCALE)
C
C       NORMALIZE INTERMEDIATE RESULT.
C
        TMP=DASUM(Z,N)
C
C       ERROR EXIT. IF LENGTH(Z)=0, SET FLAG AND RETURN.
C
        IF(TMP.EQ.0D0)THEN
          INFO=8
          RETURN
        ENDIF
C
        CALL DSCAL(Z,N,1D0/TMP)
C
C       SOLVE R*Z=Y.
C
        CALL DTRUSL(R,LDR,N,Z,S)
        IPOS(J)=IDMAX(Z,N)
   50 CONTINUE
C
C     TWO SQUARE ROOTS ARE COMPUTED BEFORE MULTIPLICATION IN ORDER TO
C     AVOID OVERFLOW.
C
      SCALE=DSQRT(DABS(SCALE)/TMP)
      TMP=DNRM2(N,Z,1)
C
C     ERROR EXIT. IF LENGTH(Z)=0, SET FLAG AND RETURN.
C
      IF(TMP.EQ.0D0)THEN
        INFO=8
        RETURN
      ENDIF
C
      DELTA=SCALE*DSQRT(DABS(S)/TMP)
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DORTHO(Q,LDQ,M,N,W,S,V,RCOND,WK,INFO)
      INTEGER LDQ,M,N,INFO
      DOUBLE PRECISION Q(LDQ,N),S(N+1),W(M),V(M),WK(M+N),RCOND
C
C     THE MATRIX Q(1:M,1:N), M>=N, IS ASSUMED TO HAVE ORTHONORMAL
C     COLUMNS. DORTHO ORTHOGONALIZES W AGAINST THE COLUMNS OF Q AND
C     RETURNS THE VECTOR SO OBTAINED IN V. IF M>N THEN V IS OF UNIT
C     LENGTH. IF M=N THEN V=0. S(1:N) IS A VECTOR OF FOURIER
C     COEFFICIENTS, I.E. S(1:N)=Q'*W WHERE THE PRIME DENOTES
C     TRANSPOSITION. S(N+1) IS THE DISTANCE FROM W TO THE RANGE OF Q.
C     IF M=N THEN S(N+1)=0.
C
C     ON ENTRY
C
C        Q       REAL*8 (LDQ,N)
C                Q CONTAINS ORTHONORMAL COLUMNS.
C
C        LDQ     INTEGER
C                LDQ IS THE LEADING DIMENSION OF THE ARRAY Q.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX Q.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF THE MATRIX Q. N MUST
C                NOT BE LARGER THAN M.
C
C        W       REAL*8 (M)
C                W IS ORTHOGONALIZED AGAINST THE COLUMNS OF Q.
C
C        RCOND   REAL*8
C                RCOND IS A BOUND FOR THE RECIPROCAL CONDITION NUMBER
C                RKAPPA OF THE M BY N+1 MATRIX OBTAINED BY APPENDING
C                THE COLUMN W/LENGHT(W) TO Q. THE RECIPROCAL CONDITION
C                NUMBER IS WITH RESPECT TO THE LEAST SQUARES NORM. THE
C                COMPUTATIONS ARE TERMINATED AND W IS NOT
C                ORTHOGONALIZED AGAINST Q, IF RCOND>RKAPPA. RCOND=0
C                PREVENTS EARLY TERMINATION.
C
C        WK      REAL*8 (M+N)
C                WK IS A SCRATCH ARRAY.
C
C     ON RETURN
C
C        V       REAL*8 (M)
C                IF W DOES NOT LIE IN THE RANGE OF Q NUMERICALLY,
C                THEN V IS A UNIT VECTOR SUCH THAT RANGE(Q,V)=
C                RANGE(Q,W), AND V IS ORTHOGONAL TO THE COLUMNS OF Q.
C                IF W LIES IN THE RANGE OF Q NUMERICALLY, THEN V=0.
C
C        S       REAL*8 (N+1)
C                S(1:N) CONTAINS THE FOURIER COEFFICIENTS OF W WITH
C                RESPECT TO THE COLUMNS OF Q, I.E. S(1:N)=Q'*W.
C                S(N+1) IS THE DISTANCE FROM W TO THE RANGE OF Q.
C
C        W       W IS THE SAME AS ON ENTRY.
C
C        RCOND   RCOND IS THE RECIPROCAL CONDITION NUMBER OF THE M BY
C                N+1 MATRIX OBTAINED BY APPENDING COLUMN W/LENGTH(W)
C                TO Q. THE CONDITION NUMBER IS WITH RESPECT TO THE
C                LEAST SQUARES NORM.
C
C        INFO    INTEGER
C
C                INFO=0    W DOES NOT LIE IN THE RANGE OF Q
C                          NUMERICALLY.
C
C                INFO=1    W WAS FOUND TO LIE IN THE RANGE OF Q
C                          NUMERICALLY. ON RETURN S(1:N)=Q'*W,
C                          S(N+1)=0, AND V=0.
C
C                INFO=2    THE RECIPROCAL CONDITION NUMBER OF THE M BY
C                          N+1 MATRIX OBTAINED BY APPENDING THE COLUMN
C                          W/LENGTH(W) TO Q IS SMALLER THAN THE ENTRY
C                          VALUE OF RCOND. ON RETURN RCOND CONTAINS
C                          THIS RECIPROCAL CONDITION NUMBER. S AND V
C                          ARE UNDEFINED.
C
C     DORTHO USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     LINPACK BLAS: DNRM2
C     OTHER BLAS: DAPX,DATPX,DAXPY,DCOPY,DSCAL,DZERO
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION DNRM2,SIGMAV,SIGMAW,SIGOLD,SIMAX,SIMIN,TMP
C
C     IF M=N, SET FLAG, S(1:N)=Q'*W, S(N+1)=0, V=0 AND RETURN.
C
      IF(M.EQ.N)THEN
        INFO=1
        CALL DATPX(Q,LDQ,M,N,W,S)
        S(N+1)=0D0
        CALL DZERO(V,M)
        RCOND=0D0
        RETURN
      ENDIF
C
C     THE LENGTH OF W IS STORED IN SIGMAW. THE LENGTH OF THE CURRENT
C     VECTOR V IS STORED IN SIGMAV. AFTER UPDATING V, THE PREVIOUS
C     VALUE OF SIGMAV IS STORED IN SIGOLD.
C
      SIGMAW=DNRM2(M,W,1)
C
C     IF W=0, SET FLAG, S=V=0 AND RETURN.
C
      IF(SIGMAW.EQ.0D0)THEN
        INFO=1
        CALL DZERO(S,N+1)
        CALL DZERO(V,M)
        RCOND=0D0
        RETURN
      ENDIF
C
C     STORE SCALED W IN V.
C
      CALL DCOPY(W,M,V)
      CALL DSCAL(V,M,1D0/SIGMAW)
C
C     ORTHOGONALIZATION OF V.
C
      CALL DATPX(Q,LDQ,M,N,V,S)
      SIMAX=DSQRT(1D0+DNRM2(N,S,1))
      CALL DAPX(Q,LDQ,M,N,S,WK(N+1))
      CALL DAXPY(WK(N+1),M,-1D0,V)
      SIGMAV=DNRM2(M,V,1)
      SIMIN=SIGMAV/SIMAX
C
C     SIMAX = LARGEST SINGULAR VALUE OF MATRIX OBTAINED BY APPENDING
C     THE COLUMN W/LENGHT(W) TO MATRIX Q.
C     SIMIM = SMALLEST SINGULAR VALUE OF MATRIX OBTAINED BY APPENDING
C     THE COLUMN W/LENGHT(W) TO MATRIX Q.
C     RETURN IF PROBLEM TOO ILL-CONDITIONED.
C
      IF(SIMIN/SIMAX.LT.RCOND)THEN
        INFO=2
        RCOND=SIMIN/SIMAX
        RETURN
      ENDIF
      RCOND=SIMIN/SIMAX
C
C     CHECK IF V IS ORTHOGONAL TO THE COLUMNS OF Q. THE TEST USED
C     IS DISCUSSED IN PARLETT, THE SYMMETRIC EIGENVALUE PROBLEM,
C     PRENTICE-HALL, ENGLEWOOD CLIFFS, N.J., 1980, P. 107.
C
      IF(SIGMAV.GT.0.717D0)THEN
        CALL DSCAL(V,M,1D0/SIGMAV)
        CALL DSCAL(S,N,SIGMAW)
        S(N+1)=SIGMAV*SIGMAW
        INFO=0
        RETURN
      ENDIF
C
C     ONE REORTHOGONALIZATION.
C
      SIGOLD=SIGMAV
      CALL DATPX(Q,LDQ,M,N,V,WK)
      CALL DAXPY(WK,N,1D0,S)
      CALL DAPX(Q,LDQ,M,N,WK,WK(N+1))
      CALL DAXPY(WK(N+1),M,-1D0,V)
      SIGMAV=DNRM2(M,V,1)
C
C     CHECK IF V IS ORTHOGONAL AGAINST THE COLUMNS OF Q.
C
      IF(SIGMAV.GT.0.717D0*SIGOLD)THEN
        CALL DSCAL(V,M,1D0/SIGMAV)
        CALL DSCAL(S,N,SIGMAW)
        S(N+1)=SIGMAV*SIGMAW
        INFO=0
        RETURN
      ENDIF
C
C     W LIES IN SPAN(Q) NUMERICALLY.
C
      INFO=1
      CALL DSCAL(S,N,SIGMAW)
      S(N+1)=0D0
      CALL DZERO(V,M)
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DORTHX(Q,LDQ,M,N,J,RHO,V,WK,INFO)
      INTEGER LDQ,M,N,INFO,J
      DOUBLE PRECISION Q(LDQ,N),V(M),WK(M+N),RHO
C
C     THE MATRIX Q(1:M,1:N), M>N, IS ASSUMED TO HAVE ORTHONORMAL
C     COLUMNS. DORTHX ORTHOGONALIZES THE JTH AXIS VECTOR AGAINST THE
C     COLUMNS OF Q AND RETURNS A UNIT VECTOR SO OBTAINED IN V. RHO IS
C     IS THE DISTANCE FROM THE JTH AXIS VECTOR TO THE RANGE OF Q.
C     DORTHX IS A SIMPLIFIED AND FASTER VERSION OF SUBROUTINE DORTHO.
C
C     ON ENTRY
C
C        Q       REAL*8 (LDQ,N)
C                Q CONTAINS ORTHONORMAL COLUMNS.
C
C        LDQ     INTEGER
C                LDQ IS THE LEADING DIMENSION OF THE ARRAY Q.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX Q.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF THE MATRIX Q. N MUST
C                NOT BE LARGER THAN M.
C
C        J       INTEGER
C                THE JTH AXIS VECTOR IS ORTHOGONALIZED AGAINST THE
C                COLUMNS OF Q.
C
C        WK      REAL*8 (M+N)
C                WK IS A SCRATCH ARRAY.
C
C     ON RETURN
C
C        V       REAL*8 (M)
C                IF THE JTH AXIS VECTOR DOES NOT LIE IN THE RANGE OF Q
C                NUMERICALLY, THEN V IS A UNIT VECTOR SUCH THAT
C                RANGE(Q,V)=RANGE(Q,(JTH AXIS VECTOR)), AND V IS
C                ORTHOGONAL TO THE COLUMNS OF Q. IF THE JTH AXIS
C                VECTOR LIES IN THE RANGE OF Q NUMERICALLY, THEN V IS
C                UNDEFINED.
C
C        RHO     REAL*8
C                RHO IS THE DISTANCE FROM THE JTH AXIS VECTOR TO THE
C                RANGE OF Q.
C
C        INFO    INTEGER
C
C                INFO=0    SUCCESSFUL TERMINATION. THE JTH AXIS VECTOR
C                          DOES NOT LIE IN THE RANGE OF Q NUMERICALLY.
C
C                INFO=8    ERROR EXIT. THE JTH AXIS VECTOR WAS FOUND
C                          TO LIE IN THE RANGE OF Q NUMERICALLY. ON
C                          RETURN V IS UNDEFINED.
C
C     DORTHX USES THE FOLLOWING FUNCTIONS AND SUBROUTINES
C
C     LINPACK BLAS: DNRM2
C     OTHER BLAS: DAPX,DATPX,DAXPY,DCOPY,DSCAL,DZERO
C
C     INTERNAL VARIABLES
C
      INTEGER I
      DOUBLE PRECISION DNRM2,RHOOLD,TMP
C
      CALL DZERO(V,M)
      V(J)=1D0
C
C     ORTHOGONALIZATION OF V. LOOP 10 IS ELIGIBLE FOR VECTORIZATION,
C     BUT DOES NOT VECTORIZE WITHOUT SPECIAL COMPILER DIRECTIVES.
C
      DO 10 I=1,N
        WK(I)=-Q(J,I)
10    CONTINUE
      CALL DAPX(Q,LDQ,M,N,WK,V)
      V(J)=V(J)+1D0
      RHO=DNRM2(M,V,1)
C
C     THE LENGTH OF THE CURRENT VECTOR V IS STORED IN RHO. AFTER
C     UPDATING V, THE PREVIOUS VALUE OF RHO IS STORED IN RHOOLD.
C
C     CHECK IF V IS ORTHOGONAL AGAINST THE COLUMNS OF Q.
C
      IF(RHO.GT.0.717D0)THEN
        CALL DSCAL(V,M,1D0/RHO)
        INFO=0
        RETURN
      ENDIF
C
C     ONE REORTHOGONALIZATION.
C
      RHOOLD=RHO
      CALL DATPX(Q,LDQ,M,N,V,WK)
      CALL DAPX(Q,LDQ,M,N,WK,WK(N+1))
      CALL DAXPY(WK(N+1),M,-1D0,V)
      RHO=DNRM2(M,V,1)
C
C     CHECK IF V IS ORTHOGONAL AGAINST THE COLUMNS OF Q.
C
      IF(RHO.GT.0.717D0*RHOOLD)THEN
        CALL DSCAL(V,M,1D0/RHO)
        INFO=0
        RETURN
      ENDIF
C
C     ERROR EXIT. THE JTH AXIS VECTOR LIES IN THE RANGE OF Q
C     NUMERICALLY.
C
      RHO=0D0
      INFO=8
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DTRLSL(T,LDT,N,B,SCALE)
C
C     DTRLSL SOLVES THE LOWER TRIANGULAR NONSINGULAR LINEAR SYSTEM
C     OF EQUATIONS WITH N BY N MATRIX TRANSPOSE(T) AND RIGHT HAND SIDE
C     B. THE MATRIX IS ASSUMED TO BE ILL-CONDITIONED, AND FREQUENT
C     RESCALINGS ARE CARRIED OUT IN ORDER TO AVOID OVERFLOW.
C
C     ON ENTRY
C
C        T       REAL*8 (LDA,N)
C                T IS AN N BY N UPPER TRIANGULAR NONSINGULAR MATRIX.
C
C        LDT     INTEGER
C                LDT IS THE LEADING DIMENSION OF THE ARRAY T.
C
C        N       INTEGER
C                N IS THE NUMBER OF ROWS AND COLUMNS OF THE MATRIX T.
C
C        B       REAL*8 (N)
C                B IS THE RIGHT HAND SIDE VECTOR.
C
C     ON RETURN
C
C        B       B CONTAINS THE SCALED SOLUTION.
C
C        SCALE   REAL*8
C                SCALE IS A SCALING FACTOR INTRODUCED IN ORDER TO
C                AVOID OVERFLOW. THE SOLUTION OF THE GIVEN SYSTEM
C                OF EQUATIONS IS B/SCALE. SCALE MAY BE NEGATIVE.
C
C     DTRLSL USES THE DOUBLE PRECISION FUNCTION DDOT.
C
      INTEGER LDT,N,J
      DOUBLE PRECISION T(LDT,N),B(N),SCALE,TMP,S,DDOT
C
      SCALE=1D0
      IF(DABS(B(1)).GT.DABS(T(1,1)))THEN
        S=T(1,1)/B(1)
        CALL DSCAL(B(2),N-1,S)
        B(1)=1D0
        SCALE=S*SCALE
      ELSE
        B(1)=B(1)/T(1,1)
      ENDIF
      IF(N.EQ.1)RETURN
C
      DO 10 J=2,N
        TMP=B(J)-DDOT(T(1,J),B,J-1)
        IF(DABS(TMP).GT.DABS(T(J,J)))THEN
          S=T(J,J)/TMP
          CALL DSCAL(B,N,S)
          B(J)=1D0
          SCALE=S*SCALE
        ELSE
          B(J)=TMP/T(J,J)
        ENDIF
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DTRUSL(U,LDU,N,B,SCALE)
C
C     DTRUSL SOLVES THE UPPER TRIANGULAR NONSINGULAR LINEAR SYSTEM
C     OF EQUATIONS WITH N BY N MATRIX U AND RIGHT HAND SIDE B. THE
C     MATRIX IS ASSUMED TO BE ILL-CONDITIONED, AND FREQUENT RESCALINGS
C     ARE CARRIED OUT IN ORDER TO AVOID OVERFLOW.
C
C     ON ENTRY
C
C        U       REAL*8 (LDA,N)
C                U IS AN N BY N UPPER TRIANGULAR NONSINGULAR MATRIX.
C
C        LDU     INTEGER
C                LDU IS THE LEADING DIMENSION OF THE ARRAY U.
C
C        N       INTEGER
C                N IS THE NUMBER OF ROWS AND COLUMNS OF THE MATRIX U.
C
C        B       REAL*8 (N)
C                B IS THE RIGHT HAND SIDE VECTOR.
C
C     ON RETURN
C
C        B       B CONTAINS THE SCALED SOLUTION.
C
C        SCALE   REAL*8
C                SCALE IS A SCALING FACTOR INTRODUCED IN ORDER TO
C                AVOID OVERFLOW. THE SOLUTION OF THE GIVEN SYSTEM
C                OF EQUATIONS IS B/SCALE. SCALE MAY BE NEGATIVE.
C
C     DTRUSL USES THE SUBROUTINE DAXPY.
C
      INTEGER LDU,N,J
      DOUBLE PRECISION U(LDU,N),B(N),SCALE,TMP
C
      SCALE=1D0
      IF(DABS(B(N)).GT.DABS(U(N,N)))THEN
        TMP=U(N,N)/B(N)
        CALL DSCAL(B,N-1,TMP)
        B(N)=1D0
        SCALE=TMP*SCALE
      ELSE
        B(N)=B(N)/U(N,N)
      ENDIF
      IF(N.EQ.1)RETURN
C
      DO 10 J=N-1,1,-1
        TMP=-B(J+1)
        CALL DAXPY(U(1,J+1),J,TMP,B)
        IF(DABS(B(J)).GT.DABS(U(J,J)))THEN
          TMP=U(J,J)/B(J)
          CALL DSCAL(B,N,TMP)
          B(J)=1D0
          SCALE=TMP*SCALE
        ELSE
          B(J)=B(J)/U(J,J)
        ENDIF
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
C BLAS THAT VECTORIZE WELL ON THE IBM 3090-200VF COMPUTER WITHOUT
C SPECIAL COMPILER DIRECTIVES.
C---------------------------------------------------------------------
      SUBROUTINE DADFGR(X,Y,C,S)
C
C     DADFGR DEFINES AND APPLIES THE SYMMETRIC GIVENS REFLECTOR
C
C
C             ( C   S )
C         G = (       )            C*C+S*S=1,
C             ( S  -C ),
C
C
C     FOR WHICH ( X, Y ) * G = ( Z, 0 ). REPLACES ( X, Y ) BY ( Z, 0 ).
C
      DOUBLE PRECISION C,MU,S,T,X,X1,Y,Y1
C
      IF(Y.EQ.0D0)THEN
        C=1D0
        S=0D0
        RETURN
      ENDIF
C
      MU=DMAX1(DABS(X),DABS(Y))
      X1=X/MU
      Y1=Y/MU
      T=DSIGN(MU*DSQRT(X1*X1+Y1*Y1),X)
      C=X/T
      S=Y/T
      X=T
      Y=0D0
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DAPLRC(C,S,N,X,Y)
C
C     DAPLRC APPLIES A GIVENS REFLECTOR G DEFINED BY C AND S, WHERE
C     G HAS BEEN DETERMINED BY SUBROUTINE DADFGR. WHEN CALLED DAPLRC
C     REPLACES THE TWO COLUMN MATRIX ( X, Y )  BY ( X, Y ) * G.
C
      INTEGER I,N
      DOUBLE PRECISION C,S,TMP,X(N),Y(N)
C
C     LOOP VECTORIZES.
C
      DO 10 I=1,N
        TMP=C*X(I)+S*Y(I)
        Y(I)=S*X(I)-C*Y(I)
        X(I)=TMP
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DAPLRR(C,S,A,LDA,N,IRW,L)
C
C     DAPLRR APPLIES A GIVENS REFLECTOR G DEFINED BY C AND S, WHERE
C     G HAS BEEN DETERMINED BY SUBROUTINE DADFGR. WHEN CALLED DAPLRR
C     REPLACES THE TWO ROW MATRIX ( A(IRW,J), A(IRW+1,J) )', L<=J<=N,
C     BY G * ( A(IRW,J), A(IRW+1,J) )', L<=J<=N.
C
C     ON ENTRY
C
C        C, S    REAL*8
C                C AND S DEFINE A GIVENS REFLECTOR G.
C
C        A       REAL*8 (LDA,N)
C                G IS APPLIED TO (PARTS OF) TWO ROWS OF THE MATRIX A.
C
C        LDA     INTEGER
C                LDA IS THE LEADING DIMENSION OF THE ARRAY A.
C
C        N       INTEGER
C                N IS THE NUMBER OF ROWS AS WELL AS THE NUMBER OF
C                COLUMNS OF THE MATRIX A.
C
C        IRW     INTEGER
C                IRW AND IRW+1 ARE THE ROW INDICES OF THE ROWS
C                OF THE MATRIX A TO WHICH G IS APPLIED.
C
C        L       INTEGER
C                G IS APPLIED TO THE 2 BY N-L+1 MATRIX
C                ( A(IRW,J), A(IRW+1,J) ), J=L,L+1,...,N.
C
C     ON RETURN
C
C        A       ( A(IRW,J), A(IRW+1,J) ), J=L,L+1,...,N, IS OBTAINED
C                BY MULTIPLYING ROWS IRW AND IRW+1, AND COLUMNS
C                L,L+1,...,N OF THE MATRIX A GIVEN ON ENTRY BY THE
C                MATRIX G.
C
      INTEGER I,IRW,L,LDA,N
      DOUBLE PRECISION C,S,TMP,A(LDA,N)
C
C     LOOP IS ELIGIBLE FOR VECTORIZATION, BUT DOES NOT VECTORIZE
C     WITHOUT SPECIAL COMPILER DIRECTIVES.
C
      DO 10 I=L,N
        TMP=C*A(IRW,I)+S*A(IRW+1,I)
        A(IRW+1,I)=S*A(IRW,I)-C*A(IRW+1,I)
        A(IRW,I)=TMP
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DAPLRV(C,S,U,USTRID,V,VSTRID,L)
C     DAPLRV APPLIES A GIVENS REFLECTOR G DEFINED BY C AND S, WHERE
C     G HAS BEEN DETERMINED BY SUBROUTINE DADFGR. WHEN CALLED DAPLRR
C     REPLACES THE TWO VECTORS ( U(1+J*USTRID), V(1+J*VSTRID) ) BY
C     G * ( U(1+J*USTRID), V(1+J*VSTRID) ), J=0,1,...L-1. DAPLRV IS
C     USED TO APPLY GIVENS REFLECTORS TO NON-ADJACENT ROWS OF A MATRIX.
C
C     ON ENTRY
C
C        C, S    REAL*8
C                C AND S DEFINE A GIVENS REFLECTOR G.
C
C        U       REAL*8 (1+(L-1)*USTRID)
C                G IS APPLIED TO U.
C
C        USTRID  INTEGER
C                USTRID IS THE STRIDE OF U.
C
C        V       REAL*8 (1+(L-1)*VSTRID)
C                G BE APPLIED TO V.
C
C        VSTRID  INTEGER
C                VSTRID IS THE STRIDE OF V.
C
C        L       INTEGER
C                L IS THE NUMBER OF TIMES G IS TO BE APPLIED.
C
C     ON RETURN
C
C        U,V     U(1+J*USTRID), V(1+J*VSTRID), J=0,1,...,L-1, ARE
C                COMPUTED BY APPLYING G TO THE CORRESPONDING ELEMENTS
C                OF U AND V GIVEN ON ENTRY.
C
      INTEGER I,L,L1,USTRID,VSTRID
      DOUBLE PRECISION C,S,TMP,U(1),V(1)
C
C     LOOP IS CONSIDERED RECURSIVE AND NOT VECTORIZABLE.
C
      L1=L-1
      DO 10 I=0,L1
        TMP=C*U(1+I*USTRID)+S*V(1+I*VSTRID)
        V(1+I*VSTRID)=S*U(1+I*USTRID)-C*V(1+I*VSTRID)
        U(1+I*USTRID)=TMP
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DAPX(A,LDA,M,N,X,Y)
C
C     DAPX COMPUTES Y:=A*X, WHERE A IS A MATRIX AND X IS A VECTOR.
C
C     ON ENTRY
C
C        A       REAL*8 (LDA,N)
C                A IS THE MATRIX THAT MULTIPLIES THE VECTOR X.
C
C        LDA     INTEGER
C                LDA IS THE LEADING DIMENSION OF THE ARRAY A.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX A.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF THE MATRIX A.
C
C        X       REAL*8 (N)
C                X IS MULTIPLIED BY THE MATRIX A.
C
C     ON RETURN
C
C        Y       REAL*8 (M)
C                Y CONTAINS THE RESULT OF THE MULTIPLICATION OF THE
C                MATRIX A TIMES THE VECTOR X.
C
      INTEGER LDA,M,N,I,J
      DOUBLE PRECISION A(LDA,N),X(N),Y(M),ACC
C
C     OUTER LOOP VECTORIZES.
C
      DO 10 I=1,M
        ACC=0D0
        DO 20 J=1,N
          ACC=ACC+A(I,J)*X(J)
   20   CONTINUE
        Y(I)=ACC
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION DASUM(X,N)
C
C     DASUM COMPUTES THE SUM OF THE ABSOLUTE VALUES OF THE COMPONENTS
C     OF X(1:N).
C
      INTEGER I,N
      DOUBLE PRECISION X(N),ACC
C
C     LOOP VECTORIZES.
C
      ACC=0D0
      DO 10 I=1,N
        ACC=ACC+DABS(X(I))
   10 CONTINUE
      DASUM=ACC
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DATPX(A,LDA,M,N,X,Y)
C
C     DATPX COMPUTES Y:=A'*X, WHERE A' DENOTES THE TRANSPOSE OF THE
C     MATRIX A AND X IS A VECTOR.
C
C     ON ENTRY
C
C        A       REAL*8 (LDA,N)
C                A IS THE MATRIX WHOSE TRANSPOSE MULTIPLIES THE VECTOR
C                X.
C
C        LDA     INTEGER
C                LDA IS THE LEADING DIMENSION OF THE ARRAY A.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX A.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF THE MATRIX A.
C
C        X       REAL*8 (M)
C                X IS MULTIPLIED BY THE TRANSPOSE OF THE MATRIX A.
C
C     ON RETURN
C
C        Y       REAL*8 (N)
C                Y CONTAINS THE RESULT OF THE MULTIPLICATION OF THE
C                TRANSPOSE OF MATRIX A TIMES THE VECTOR X.
C
C     DATPX USES THE DOUBLE PRECISION FUNCTION DDOT.
C
      INTEGER LDA,M,N,I
      DOUBLE PRECISION A(LDA,N),X(M),Y(N),DDOT
C
      DO 10 I=1,N
        Y(I)=DDOT(A(1,I),X,M)
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DAXPY(X,N,ALPHA,Y)
C
C     DAXPY COMPUTES Y:=Y+ALPHA*X, WHERE X AND Y ARE VECTORS, AND
C     ALPHA IS A CONSTANT.
C
      INTEGER N,I
      DOUBLE PRECISION X(N),Y(N),ALPHA
C
C     LOOP VECTORIZES.
C
      DO 10 I=1,N
        Y(I)=Y(I)+ALPHA*X(I)
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DCOPY(X,N,Y)
C
C     DCOPY COPIES A VECTOR X TO A VECTOR Y.
C
      INTEGER I,N
      DOUBLE PRECISION X(N),Y(N)
C
C     LOOP VECTORIZES.
C
      DO 10 I=1,N
        Y(I)=X(I)
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION DDOT(X,Y,N)
C
C     DDOT COMPUTES THE INNER PRODUCT OF X(1:N) AND Y(1:N).
C
      INTEGER N,I
      DOUBLE PRECISION X(N),Y(N),ACC
C
C     LOOP VECTORIZES.
C
      ACC=0D0
      DO 10 I=1,N
        ACC=ACC+X(I)*Y(I)
   10 CONTINUE
      DDOT=ACC
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DMINRW(A,LDA,M,N,X,I)
C
C     DMINRW DETERMINES A ROW OF THE MATRIX A WITH LEAST LEAST SQUARES
C     NORM.
C
C     ON ENTRY
C
C        A       REAL*8 (LDA,N)
C                THE MATRIX WHOSE ROW WITH LEAST LEAST SQUARES NORM IS
C                DETERMINED.
C
C        LDA     INTEGER
C                LDA IS THE LEADING DIMENSION OF THE ARRAY A.
C
C        M       INTEGER
C                M IS THE NUMBER OF ROWS OF THE MATRIX A.
C
C        N       INTEGER
C                N IS THE NUMBER OF COLUMNS OF THE MATRIX A.
C
C     ON RETURN
C
C        X       REAL*8 (M)
C                X CONTAINS THE LEAST SQUARES NORM OF THE ROWS OF THE
C                MATRIX A.
C
C        I       INTEGER
C                I CONTAINS THE ROW NUMBER OF THE FIRST ROW OF THE
C                MATRIX A OF MINIMAL LEAST SQUARES NORM.
C
      INTEGER LDA,M,N,I,J,K
      DOUBLE PRECISION A(LDA,N),X(M),ACC,T
C
C     OUTER LOOP VECTORIZES.
C
      DO 10 J=1,M
        ACC=0D0
        DO 20 K=1,N
          ACC=ACC+A(J,K)*A(J,K)
   20   CONTINUE
        X(J)=ACC
   10 CONTINUE
C
C     DETERMINE SMALLEST COMPONENT OF X. LOOP DOES NOT VECTORIZE.
C
      I=1
      T=X(1)
      DO 30 J=2,M
        IF(X(J).GE.T)GOTO 30
        T=X(J)
        I=J
   30 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DSCAL(X,N,ALPHA)
C
C     DSCAL COMPUTES A CONSTANT TIMES A VECTOR.
C
      INTEGER N,I
      DOUBLE PRECISION X(N),ALPHA
C
C     LOOP VECTORIZES.
C
      DO 10 I=1,N
        X(I)=ALPHA*X(I)
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DSHFTD(X,N)
C
C     DSHFTD SHIFTS A VECTOR X(2:N) DOWNWARDS 1 POSITION.
C
      INTEGER I,N
      DOUBLE PRECISION X(N)
C
C     LOOP VECTORIZES.
C
      DO 10 I=2,N
        X(I-1)=X(I)
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DSHFTU(X,N)
C
C     DSHFTU SHIFTS A VECTOR X(1:N) UPWARDS 1 POSITION.
C
      INTEGER I,N
      DOUBLE PRECISION X(N+1)
C
C     LOOP VECTORIZES.
C
      DO 10 I=N,1,-1
        X(I+1)=X(I)
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      SUBROUTINE DZERO(X,N)
C
C     DZERO SETS A VECTOR TO ZERO.
C
      INTEGER I,N
      DOUBLE PRECISION X(N)
C
C     LOOP VECTORIZES.
C
      DO 10 I=1,N
        X(I)=0D0
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
      INTEGER FUNCTION IDMAX(X,N)
C
C     IDMAX DETERMINES THE INDEX OF AN ELEMENT OF VECTOR X(1:N) OF
C     LARGEST ABSOLUTE VALUE. IF THERE ARE SEVERAL SUCH ELEMENTS, THE
C     SMALLEST POSSIBLE INDEX IS RETURNED.
C
      INTEGER I,N
      DOUBLE PRECISION X(N),TMP
C
C     LOOP CONSIDERED RECURSIVE AND DOES NOT VECTORIZE.
C
      IF(N.LE.0)RETURN
      IDMAX=1
      IF(N.EQ.1)RETURN
      TMP=DABS(X(1))
      DO 10 I=2,N
        IF(DABS(X(I)).LE.TMP)GOTO 10
        IDMAX=I
        TMP=DABS(X(I))
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------
C LINPACK BLAS DNRM2.
C---------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
      INTEGER          NEXT
      DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
      DATA   ZERO, ONE /0.0D0, 1.0D0/
C
C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
C     INCREMENT INCX .
C     IF    N .LE. 0 RETURN WITH RESULT = 0.
C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C           C.L.LAWSON, 1978 JAN 08
C
C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
C     HOPEFULLY APPLICABLE TO ALL MACHINES.
C         CUTLO = MAXIMUM OF  DSQRT(U/EPS)  OVER ALL KNOWN MACHINES.
C         CUTHI = MINIMUM OF  DSQRT(V)      OVER ALL KNOWN MACHINES.
C     WHERE
C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
C
C     BRIEF OUTLINE OF ALGORITHM..
C
C     PHASE 1    SCANS ZERO COMPONENTS.
C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C     VALUES FOR CUTLO AND CUTHI..
C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
C                   UNIVAC AND DEC AT 2**(-103)
C                   THUS CUTLO = 2**(-51) = 4.44089E-16
C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C                   THUS CUTHI = 2**(63.5) = 1.30438E19
C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
      DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C
      IF(N .GT. 0) GO TO 10
         DNRM2  = ZERO
         GO TO 300
C
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
      NN = N * INCX
C                                                 BEGIN MAIN LOOP
      I = 1
   20    GO TO NEXT,(30, 50, 70, 110)
   30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
      ASSIGN 50 TO NEXT
      XMAX = ZERO
C
C                        PHASE 1.  SUM IS ZERO
C
   50 IF( DX(I) .EQ. ZERO) GO TO 200
      IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
C
C                                PREPARE FOR PHASE 2.
      ASSIGN 70 TO NEXT
      GO TO 105
C
C                                PREPARE FOR PHASE 4.
C
  100 I = J
      ASSIGN 110 TO NEXT
      SUM = (SUM / DX(I)) / DX(I)
  105 XMAX = DABS(DX(I))
      GO TO 115
C
C                   PHASE 2.  SUM IS SMALL.
C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
   70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
C
C                     COMMON CODE FOR PHASES 2 AND 4.
C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
C
  110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
         SUM = ONE + SUM * (XMAX / DX(I))**2
         XMAX = DABS(DX(I))
         GO TO 200
C
  115 SUM = SUM + (DX(I)/XMAX)**2
      GO TO 200
C
C
C                  PREPARE FOR PHASE 3.
C
   75 SUM = (SUM * XMAX) * XMAX
C
C
C     FOR REAL OR D.P. SET HITEST = CUTHI/N
C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
C
   85 HITEST = CUTHI/FLOAT( N )
C
C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
C
      DO 95 J =I,NN,INCX
      IF(DABS(DX(J)) .GE. HITEST) GO TO 100
   95    SUM = SUM + DX(J)**2
      DNRM2 = DSQRT( SUM )
      GO TO 300
C
  200 CONTINUE
      I = I + INCX
      IF ( I .LE. NN ) GO TO 20
C
C              END OF MAIN LOOP.
C
C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
      DNRM2 = XMAX * DSQRT(SUM)
  300 CONTINUE
      RETURN
      END

		
.

NEW PAGES:

[ODDNUGGET]

[GOPHER]