SUBROUTINE CDSVDC ( * * inputs * : M, N, MP, NP, RV1, * * input/output * : A, * * outputs * : W, V, STATUS) * * Module number: * * Module name: * * Keyphrase: * ---------- * perform the sigular value decomposition * * Description: * ------------ * adapted from Numerical Recipes by William Press et. al., p. 60, * subroutine SVDCMP * Double precision version * given a matrix a, with logical dimensions M rows by N columns and physical * dimensions MP rows by NP columns, this routine computes its singular * decomposition, A = U * W * tr(V). The matrix U replaces A on output, the * diagonal matrix of singular values W is output as a vector W. The matrix V * (not the transpose) is output as V. M must be greater or equal to N; if it * is smaller, then A should be filled up to square with rows of zeros. * * FORTRAN name: CDSVDC.FOR * * Keywords of accessed files and tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: * None * SDAS: * UMSPUT * Others: * None * * History: * -------- * Version Date Author Description * 1 12-11-88 J.-C. HSU coding *------------------------------------------------------------------------------- * *== input: * --logical dimensions of input/output * --matrix A (U) INTEGER M, N, * --physical dimensions of input/output * --matrix A (U) : MP, NP * --working area, MUST >= N DOUBLE PRECISION RV1(N) * *== input/output: * --matrix to be SV decomposed DOUBLE PRECISION A(MP, NP) * *== output: * --component matrices DOUBLE PRECISION V(NP, NP), W(NP) * --error status INTEGER STATUS * *== local: * DOUBLE PRECISION ANORM, C, F, G, H, S, SCALE, X, Y, Z CHARACTER*3 CHAR3 * --error message CHARACTER*130 CONTXT, MESS * --loop indices INTEGER ITS, ITER, I, J, K, L, NM, : STATOK PARAMETER (ITER = 30) * *=========================begin hsp.inc========================================= * --status return code INTEGER OK, ERRNUM(20) INTEGER DEST, PRIO DATA OK /0/ DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710, : 711, 712, 713, 714, 715, 716, 717, 718, 719, 720/ * --message destination and priority DATA DEST, PRIO /1, 0/ *=========================end hsp.inc=========================================== *------------------------------------------------------------------------------ * * make sure the number of rows is not less than the number of of columns * IF (M .LT. N) THEN STATUS = ERRNUM(1) CONTXT = 'number of rows in input matrix is less than ' // : 'number of columns, need to augment with rows ' : // 'of zeros' GO TO 999 END IF * * Householder reduction to bidiagonal form * G = 0.D0 SCALE = 0.D0 ANORM = 0.D0 * DO 250 I = 1, N L = I + 1 RV1(I) = SCALE * G G = 0.D0 S = 0.D0 SCALE = 0.D0 * IF (I .LE. M) THEN DO 110 K = I, M SCALE = SCALE + ABS(A(K,I)) 110 CONTINUE IF (SCALE .NE. 0.D0) THEN DO 120 K = I, M A(K,I) = A(K,I) / SCALE S = S + A(K,I) * A(K,I) 120 CONTINUE F = A(I,I) G = -SIGN(SQRT(S), F) H = F * G - S A(I,I) = F - G IF (I .NE. N) THEN DO 150 J = L, N S = 0.D0 DO 130 K = I, M S = S + A(K,I) * A(K,J) 130 CONTINUE F = S / H DO 140 K = I, M A(K,J) = A(K,J) + F * A(K,I) 140 CONTINUE 150 CONTINUE ENDIF DO 160 K = I, M A(K,I) = SCALE * A(K,I) 160 CONTINUE ENDIF ENDIF * W(I) = SCALE * G G = 0.D0 S = 0.D0 SCALE = 0.D0 * IF ((I .LE. M) .AND. (I .NE. N)) THEN DO 170 K = L, N SCALE = SCALE + ABS(A(I,K)) 170 CONTINUE IF (SCALE .NE. 0.D0) THEN DO 180 K = L, N A(I,K) = A(I,K) / SCALE S = S + A(I,K) * A(I,K) 180 CONTINUE F = A(I,L) G = -SIGN(SQRT(S), F) H = F * G - S A(I,L) = F - G DO 190 K = L, N RV1(K) = A(I,K) / H 190 CONTINUE IF (I .NE. M) THEN DO 230 J = L, M S = 0.D0 DO 210 K = L, N S = S + A(J,K) * A(I,K) 210 CONTINUE DO 220 K = L, N A(J,K) = A(J,K) + S * RV1(K) 220 CONTINUE 230 CONTINUE ENDIF DO 240 K = L, N A(I,K) = SCALE * A(I,K) 240 CONTINUE ENDIF ENDIF ANORM = MAX(ANORM, (ABS(W(I))+ABS(RV1(I)))) 250 CONTINUE * * accumulation of right-hand transformation * DO 320 I = N, 1, -1 IF (I .LT. N) THEN IF (G .NE. 0.D0) THEN DO 260 J = L, N * * double division to avoid possible underflow * V(J,I) = (A(I,J) / A(I,L)) / G 260 CONTINUE DO 290 J = L, N S = 0.D0 DO 270 K = L, N S = S + A(I,K) * V(K,J) 270 CONTINUE DO 280 K = L, N V(K,J) = V(K,J) + S * V(K,I) 280 CONTINUE 290 CONTINUE ENDIF DO 310 J = L, N V(I,J) = 0.D0 V(J,I) = 0.D0 310 CONTINUE ENDIF V(I,I) = 1.D0 G = RV1(I) L = I 320 CONTINUE * * accumulation of left-hand transformation * DO 390 I = N, 1, -1 L = I + 1 G = W(I) IF (I .LT. N) THEN DO 330 J = L, N A(I,J) = 0.D0 330 CONTINUE ENDIF IF (G .NE. 0.D0) THEN G = 1.D0 / G IF (I .NE. N) THEN DO 360 J = L, N S = 0.D0 DO 340 K = L, M S = S + A(K,I) * A(K,J) 340 CONTINUE F = (S / A(I,I)) * G DO 350 K = I, M A(K,J) = A(K,J) + F * A(K,I) 350 CONTINUE 360 CONTINUE ENDIF DO 370 J = I, M A(J,I) = A(J,I) * G 370 CONTINUE ELSE DO 380 J = I, M A(J,I) = 0.D0 380 CONTINUE ENDIF A(I,I) = A(I,I) + 1.D0 390 CONTINUE * * diagonalization of the bidiagonal form * * loop over singular values * DO 490 K = N, 1, -1 * * loop over allowed iterations * DO 480 ITS = 1, ITER * * test for splitting * DO 410 L = K, 1, -1 * * note that RV1(1) is always zero * NM = L - 1 IF ((ABS(RV1(L))+ANORM) .EQ. ANORM) GO TO 20 IF ((ABS(W(NM))+ANORM) .EQ. ANORM) GO TO 10 410 CONTINUE * * cancellation of RV(1), if L > 1 * 10 C = 0.D0 S = 1.D0 DO 430 I = L, K F = S * RV1(I) IF ((ABS(F)+ANORM) .NE. ANORM) THEN G = W(I) * * fancy way to do SQRT(F**2+G**2) to avoid overflow * H = ABS(F+G) * SQRT(1.D0 - 2.D0 * (F / (F+G)) : * (G / (F+G))) W(I) = H H = 1.D0 / H C = G * H S = -F * H DO 420 J = 1, M Y = A(J,NM) Z = A(J,I) A(J,NM) = (Y*C) + (Z*S) A(J,I) = -(Y*S) + (Z*C) 420 CONTINUE ENDIF 430 CONTINUE * 20 Z = W(K) * * convergence * IF (L .EQ. K) THEN * * singular value is made nonnegative * IF (Z .LT. 0.D0) THEN W(K) = -Z DO 440 J = 1, N V(J,K) = -V(J,K) 440 CONTINUE ENDIF GO TO 30 ENDIF * IF (ITS .EQ. ITER) THEN WRITE (CHAR3, '(I3)') ITER STATUS = ERRNUM(2) CONTXT = 'No convergence in ' // CHAR3 // : 'iterations' GO TO 999 END IF * * shift from bottom 2-by-2 minor * X = W(L) NM = K - 1 Y = W(NM) G = RV1(NM) H = RV1(K) F = ((Y-Z)*(Y+Z) + (G-H)*(G+H)) / (2.0*H*Y) * * fancy way to do SQRT(F**2+1) to avoid overflow * G = ABS(F+1.D0) * SQRT(1.D0 - 2.D0 * F / (F+1.D0) / : (F+1.D0)) F = ((X-Z)*(X+Z) + H*((Y/(F+SIGN(G,F))) - H)) / X * * next QR transformation * C = 1.D0 S = 1.D0 DO 470 J = L, NM I = J + 1 G = RV1(I) Y = W(I) H = S * G G = C * G * * fancy way to do SQRT(F**2+H**2) to avoid overflow * Z = ABS(F+H) * SQRT(1.D0 - 2.D0 * (F / (F+H)) : * (H / (F+H))) RV1(J) = Z C = F / Z S = H / Z F = (X*C) + (G*S) G = -(X*S) + (G*C) H = Y * S Y = Y * C DO 450 NM = 1, N X = V(NM,J) Z = V(NM,I) V(NM,J) = (X*C) + (Z*S) V(NM,I) = -(X*S) + (Z*C) 450 CONTINUE * * fancy way to do SQRT(F**2+H**2) to avoid overflow * Z = ABS(F+H) * SQRT(1.D0 - 2.D0 * (F / (F+H)) : * (H / (F+H))) * * rotation can be arbitrary if Z = 0 * W(J) = Z IF (Z .NE. 0.D0) THEN Z = 1.D0 / Z C = F * Z S = H * Z ENDIF F = (C*G) + (S*Y) X = -(S*G) + (C*Y) DO 460 NM = 1, M Y = A(NM,J) Z = A(NM,I) A(NM,J) = (Y*C) + (Z*S) A(NM,I) = -(Y*S) + (Z*C) 460 CONTINUE 470 CONTINUE RV1(L) = 0.D0 RV1(K) = F W(K) = X 480 CONTINUE 30 CONTINUE 490 CONTINUE * GO TO 1000 * 999 MESS = 'CDSVDC: ' // CONTXT CALL UMSPUT (MESS, DEST, PRIO, STATOK) * 1000 RETURN END