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
.