************************************************************************ * Routines from LINPACK and BLAS for the singular value decomposition, * from the Naval Surface Warfare Center library. * Public Domain. ************************************************************************ 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,*),S(*),E(*),U(LDU,*),V(LDV,*),WORK(*) 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 rows of the matrix X. C C P INTEGER. C P is the number of columns 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 left 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 DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG 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 DDOTN,T DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,DNRM2N,SCALE,SHIFT,SL,SM,SN, $ SMM1,T1,TEST,ZTEST LOGICAL WANTU,WANTV * Avoid compiler warnings DATA L, LS / 0, 0 / 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) = DNRM2N(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 DSCALN(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 = -DDOTN(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL DAXPYN(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) = DNRM2N(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 DSCALN(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 DAXPYN(N-L,E(J),X(LP1,J),1,WORK(LP1),1) 100 CONTINUE DO 110 J = LP1, P CALL DAXPYN(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 of 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 = -DDOTN(N-L+1,U(L,L),1,U(L,J),1)/U(L,L) CALL DAXPYN(N-L+1,T,U(L,L),1,U(L,J),1) 210 CONTINUE 220 CONTINUE CALL DSCALN(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 = -DDOTN(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L) CALL DAXPYN(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 DROTGN(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 DROTN(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 DROTGN(T1,F,CS,SN) S(K) = T1 F = -SN*E(K) E(K) = CS*E(K) IF (WANTU) CALL DROTN(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 DROTGN(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 DROTN(P,V(1,K),1,V(1,K+1),1,CS,SN) CALL DROTGN(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 DROTN(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 DSCALN(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 DSWAPN(P,V(1,L),1,V(1,L+1),1) IF (WANTU .AND. L .LT. N) $ CALL DSWAPN(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 ************************************************************************ SUBROUTINE DAXPYN(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(*),DY(*),DA INTEGER I,INCX,INCY,IX,IY,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 ************************************************************************ DOUBLE PRECISION FUNCTION DNRM2N ( N, DX, INCX) INTEGER NEXT DOUBLE PRECISION DX(*), 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 * Avoid compiler warnings DATA XMAX / 0.D0 / IF(N .GT. 0) GO TO 10 DNRM2N = 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/DBLE( 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 DNRM2N = 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 DNRM2N = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END ************************************************************************ DOUBLE PRECISION FUNCTION DDOTN(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(*),DY(*),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOTN = 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 DDOTN = 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 DDOTN = DTEMP RETURN END ************************************************************************ SUBROUTINE DSCALN(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(*) 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 ************************************************************************ SUBROUTINE DROTGN(DA,DB,DC,DS) C C Designed by C.L.Lawson, JPL, 1977 Sept 08 C C C Construct the Givens transformation C C ( DC DS ) C G = ( ) , DC**2 + DS**2 = 1 , C (-DS DC ) C C Which zeros the second entry of the 2-vector (DA,DB)**T . C C The quantity R = (+/-)DSQRT(DA**2 + DB**2) overwrites DA in C storage. The value of DB is overwritten by a value Z which c allows DC and DS to be recovered by the following algorithm: C If Z=1 set DC=0.D0 and DS=1.D0 C If DABS(Z) .LT. 1 set DC=DSQRT(1-Z**2) and DS=Z C If DABS(Z) .GT. 1 set DC=1/Z AND DS=DSQRT(1-DC**2) C C Normally, the subprogram DROT(N,DX,INCX,DY,INCY,DC,DS) will C next be called to apply the transformation to a 2 by N matrix. C C ------------------------------------------------------------------ C DOUBLE PRECISION DA, DB, DC, DS, U, V, R IF (DABS(DA) .LE. DABS(DB)) GO TO 10 C C *** Here DABS(DA) .GT. DABS(DB) *** C U = DA + DA V = DB / U C C Note that U and R have the sign of DA C R = DSQRT(.25D0 + V**2) * U C C Note that DC is positive C DC = DA / R DS = V * (DC + DC) DB = DS DA = R RETURN C C *** Here DABS(DA) .LE. DABS(DB) *** C 10 IF (DB .EQ. 0.D0) GO TO 20 U = DB + DB V = DA / U C C Note that U and R have the sign of DB C (R is immediately stored in DA) C DA = DSQRT(.25D0 + V**2) * U C C Note that DS is positive C DS = DB / DA DC = V * (DS + DS) IF (DC .EQ. 0.D0) GO TO 15 DB = 1.D0 / DC RETURN 15 DB = 1.D0 RETURN C C *** Here DA = DB = 0.D0 *** C 20 DC = 1.D0 DS = 0.D0 RETURN C END ************************************************************************ SUBROUTINE DROTN (N,DX,INCX,DY,INCY,C,S) C C Applies a plane rotation. C Jack Dongarra, LINPACK, 3/11/78. C DOUBLE PRECISION DX(*),DY(*),DTEMP,C,S INTEGER I,INCX,INCY,IX,IY,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 = C*DX(IX) + S*DY(IY) DY(IY) = C*DY(IY) - S*DX(IX) DX(IX) = DTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C Code for both increments equal to 1 C 20 DO 30 I = 1,N DTEMP = C*DX(I) + S*DY(I) DY(I) = C*DY(I) - S*DX(I) DX(I) = DTEMP 30 CONTINUE RETURN END ************************************************************************ SUBROUTINE DSWAPN (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(*),DY(*),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 *+++++++++++++++++++++++ End of file dsvdc.f +++++++++++++++++++++++++++