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 MN OR MM OR M<=N OR LDQM OR M<=N OR LDQ0 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 KN OR MM OR LDQM OR M<=N OR LDQN OR MN OR M0, 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