'************************************************************************ '* SOLVE A LINEAR SYSTEM AX=B BY THE SINGULAR VALUE DECOMPOSITION METHOD * '* * '* This method is recommanded when the matrix A is ill-conditionned or * '* near-singular, or there are fewer equations than unknowns. * '* --------------------------------------------------------------------- * '* SAMPLE RUN: * '* 1. Solve normal linear system of size 5 x 5: * '* Matrix A: Vector B: * '* 1 2 0 0 0 3 * '* 2 4 7 0 0 13 * '* 0 7 10 8 0 25 * '* 0 0 8 -0.75 -9 -1.75 * '* 0 0 0 -9 10 1 * '* * '* Input file contains: * '* 5 * '* 5 * '* 1.000000 2.000000 0.000000 0.000000 0.000000 * '* 2.000000 4.000000 7.000000 0.000000 0.000000 * '* 0.000000 7.000000 10.000000 8.000000 0.000000 * '* 0.000000 0.000000 8.000000 -0.750000 -9.000000 * '* 0.000000 0.000000 0.000000 -9.000000 10.000000 * '* 3.000000 * '* 13.000000 * '* 25.000000 * '* -1.750000 * '* 1.000000 * '* * '* Output file contains: * '* M = 5 * '* N = 5 * '* Matrix A: * '* 1.000000 2.000000 0.000000 0.000000 0.000000 * '* 2.000000 4.000000 7.000000 0.000000 0.000000 * '* 0.000000 7.000000 10.000000 8.000000 0.000000 * '* 0.000000 0.000000 8.000000 -0.750000 -9.000000 * '* 0.000000 0.000000 0.000000 -9.000000 10.000000 * '* Vector B: * '* 3.000000 13.000000 25.000000 -1.750000 1.000000 * '* Matrix U: * '* -0.034422 0.690358 -0.717675 0.071746 -0.044894 * '* -0.311812 -0.614665 -0.549190 0.413973 0.227984 * '* -0.663543 0.222546 0.320343 0.484542 -0.415673 * '* -0.483350 0.237931 0.181439 -0.208391 0.795873 * '* 0.477150 0.198631 0.218615 0.738425 0.373911 * '* Vector W: * '* 19.116947 2.530469 0.780714 12.539891 9.156592 * '* Matrix V: * '* -0.034422 -0.690358 -0.717675 0.071746 0.044894 * '* -0.311812 0.614665 -0.549190 0.413973 -0.227984 * '* -0.663543 -0.222546 0.320343 0.484542 0.415673 * '* -0.483350 -0.237931 0.181439 -0.208391 -0.795873 * '* 0.477150 -0.198631 0.218615 0.738425 -0.373911 * '* Solution: * '* 1.000000 1.000000 1.000000 1.000000 1.000000 * '* * '* 2. Solve linear system of size 4 x 5 (m < n): * '* Matrix A: Vector B: * '* 1 2 0 0 0 3 * '* 2 4 7 0 0 13 * '* 0 7 10 8 0 25 * '* 0 0 8 -0.75 -9 -1.75 * '* * '* Input file contains: * '* 4 * '* 5 * '* 1.000000 2.000000 0.000000 0.000000 0.000000 * '* 2.000000 4.000000 7.000000 0.000000 0.000000 * '* 0.000000 7.000000 10.000000 8.000000 0.000000 * '* 0.000000 0.000000 8.000000 -0.750000 -9.000000 * '* 3.000000 * '* 13.000000 * '* 25.000000 * '* -1.750000 * '* * '* Output file contains: * '* M = 4 * '* N = 5 * '* Matrix A: * '* 1.000000 2.000000 0.000000 0.000000 0.000000 * '* 2.000000 4.000000 7.000000 0.000000 0.000000 * '* 0.000000 7.000000 10.000000 8.000000 0.000000 * '* 0.000000 0.000000 8.000000 -0.750000 -9.000000 * '* Vector B: * '* 3.000000 13.000000 25.000000 -1.750000 * '* Matrix U: * '* -0.048526 0.080513 -0.252041 0.000000 0.963140 * '* -0.420052 0.055591 -0.869634 0.000000 -0.253383 * '* -0.769221 0.500339 0.396761 -0.000000 0.023246 * '* -0.479062 -0.860284 0.150972 -0.000000 0.087286 * '* Vector W: * '* 17.705504 9.945129 4.179107 0.000000 1.645325 * '* Matrix V: * '* -0.050190 0.019275 -0.476492 0.832543 0.277376 * '* -0.404496 0.390720 -0.288408 -0.416271 0.653651 * '* -0.816982 -0.149796 -0.218240 -0.000000 -0.512322 * '* -0.327269 0.467357 0.732420 0.364238 0.073239 * '* 0.243515 0.778527 -0.325129 -0.030353 -0.477457 * '* Solution: * '* 0.375463 1.312268 1.000000 0.726765 1.022770 * '* * '* 3. Solve linear system of size 5 x 5 (near-singular): * '* Matrix A: Vector B: * '* 1 2 0 0 0 3 * '* 2 4 1e-6 0 0 6.000001 * '* 0 7 10 8 0 25 * '* 0 0 8 -0.75 -9 -1.75 * '* 0 0 0 -9 10 1 * '* * '* Input file contains: * '* 5 * '* 5 * '* 1.000000 2.000000 0.000000 0.000000 0.000000 * '* 2.000000 4.000000 0.000001 0.000000 0.000000 * '* 0.000000 7.000000 10.000000 8.000000 0.000000 * '* 0.000000 0.000000 8.000000 -0.750000 -9.000000 * '* 0.000000 0.000000 0.000000 -9.000000 10.000000 * '* 3.000000 * '* 6.000001 * '* 25.000000 * '* -1.750000 * '* 1.000000 * '* * '* Output file contains: * '* M = 5 * '* N = 5 * '* Matrix A: * '* 1.000000 2.000000 0.000000 0.000000 0.000000 * '* 2.000000 4.000000 0.000001 0.000000 0.000000 * '* 0.000000 7.000000 10.000000 8.000000 0.000000 * '* 0.000000 0.000000 8.000000 -0.750000 -9.000000 * '* 0.000000 0.000000 0.000000 -9.000000 10.000000 * '* Vector B: * '* 3.000000 6.000000 25.000000 -1.750000 1.000000 * '* Matrix U: * '* -0.029197 0.894427 -0.435067 0.093070 -0.034671 * '* -0.058394 -0.447214 -0.870134 0.186141 -0.069341 * '* -0.649209 0.000000 0.207861 0.720418 -0.127747 * '* -0.500300 0.000000 -0.091238 -0.280151 0.814181 * '* 0.569179 0.000000 0.045304 0.599335 0.561053 * '* Vector W: * '* 18.338461 0.000000 4.279165 11.548505 8.751235 * '* Matrix V: * '* -0.007961 -0.859939 -0.508355 0.040295 -0.019809 * '* -0.263732 0.429970 -0.676684 0.517265 -0.141801 * '* -0.572267 -0.174666 0.315179 0.429750 0.598313 * '* -0.542088 -0.157890 0.309309 0.050175 -0.763559 * '* 0.555908 -0.142101 0.297764 0.737300 -0.196212 * '* Solution: * '* 0.222075 1.388963 0.841992 0.857168 0.871451 * '* * '* --------------------------------------------------------------------- * '* Reference: "Numerical Recipes by W.H. Press, B. P. Flannery, * '* S.A. Teukolsky and W.T. Vetterling, Cambridge * '* University Press, 1986" [BIBLI 08]. * '* * '* Basic Release 2.0 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* --------------------------------------------------------------------- * '* Release 2.0: Added dynamic allocations, I/O files and two more * '* examples. * '************************************************************************* 'PROGRAM TEST_SVBKSB ' LIST OF SUBROUTINES: ' =================== ' 500 Computes sqr(xa*xa + xb*xb) without destructive ' underflow or overflow. ' 600 Emulation of Fortran Function Sign. ' 1000 Solves a linear system A · X = B by Singular Value ' Decomposition Method. ' 2000 Singular Value Decomposition: Matrix A is transformed ' into product U . W . Vt (must be called before 1000). '------------------------------------------------------------------------- DEFDBL A-H, O-Z DEFINT I-N 'open input file CLS PRINT INPUT " Input data file name: ", name\$ OPEN name\$ FOR INPUT AS #1 'read size of linear system INPUT #1, M INPUT #1, N F\$ = "#####.######" 'print format NMAX = 100 'max size of utility vector TMP DIM A(M, N), U(M, N), W(N), V(N, N), B(N), X(N) DIM TMP(NMAX) 'Read matrix A FOR i = 1 TO M FOR j = 1 TO N INPUT #1, A(i, j) NEXT j NEXT i 'Read vector B FOR i = 1 TO M INPUT #1, B(i) NEXT i CLOSE #1 'open input file OPEN "tsvbksb.lst" FOR OUTPUT AS #2 PRINT #2, PRINT #2, " M = "; M PRINT #2, " N = "; N PRINT #2, " Matrix A:" FOR i = 1 TO M FOR j = 1 TO N - 1 PRINT #2, USING F\$; A(i, j); NEXT j PRINT #2, USING F\$; A(i, N) NEXT i PRINT #2, " Vector B:" FOR j = 1 TO M - 1 PRINT #2, USING F\$; B(j); NEXT j PRINT #2, USING F\$; B(M) ' Copy A in U FOR i = 1 TO M FOR j = 1 TO N U(i, j) = A(i, j) NEXT j NEXT i 'call singular value decomposition subroutine GOSUB 2000 'Call SVDCMP(U,M,N,W,V) 'print intermediate results PRINT #2, " Matrix U:" FOR i = 1 TO M FOR j = 1 TO N - 1 PRINT #2, USING F\$; U(i, j); NEXT j PRINT #2, USING F\$; U(i, N) NEXT i PRINT #2, " Vector W:" FOR j = 1 TO N - 1 PRINT #2, USING F\$; W(j); NEXT j PRINT #2, USING F\$; W(N) PRINT #2, " Matrix V:" FOR i = 1 TO N FOR j = 1 TO N - 1 PRINT #2, USING F\$; V(i, j); NEXT j PRINT #2, USING F\$; V(i, N) NEXT i 'seek highest value of W"s and set near-zero 'values to exactly zero (for near-singular cases) WMAX = 0# FOR j = 1 TO N IF W(j) > WMAX THEN WMAX = W(j) NEXT j WMIN = WMAX * .000001# FOR j = 1 TO N IF W(j) < WMIN THEN W(j) = 0# NEXT j 'call solver for SVD matrix GOSUB 1000 'Call SVBKSB(U,W,V,M,N,B,X) 'print solution PRINT #2, " Solution:" FOR j = 1 TO N PRINT #2, USING F\$; X(j); NEXT j PRINT #2, CLOSE #2 PRINT PRINT " Results in file tsvbksb.lst. " PRINT END 'of main program 'Utility functions 500 'Function Pythag(xa,xb) 'Computes sqr(xa*xa + xb*xb) without destructive underflow or overflow absa = ABS(xa) absb = ABS(xb) IF absa > absb THEN Pythag = absa * SQR(1# + (absb / absa) * (absb / absa)) ELSE IF absb = 0# THEN Pythag = 0# ELSE Pythag = absb * SQR(1# + (absa / absb) * (absa / absb)) END IF END IF RETURN 600 'Function Sign(xa,xb) IF xb < 0# THEN Sign = -ABS(xa) ELSE Sign = ABS(xa) END IF RETURN 1000 'SUBROUTINE svbksb(U,W,V,M,N,B,X) '------------------------------------------------------------------------------------------- ' Solves A · X = B for a vector X, where A is specified by the arrays U, W, V as returned by ' svdcmp. m and n are the dimensions of A, and will be equal for square matrices. B(1:m) is ' the input right-hand side. X(1:n) is the output solution vector. No input quantities are ' destroyed, so the routine may be called sequentially with different B's. '------------------------------------------------------------------------------------------- FOR j = 1 TO N 'Calculate UTB. s = 0# IF W(j) <> 0# THEN 'Nonzero result only if wj is nonzero. FOR i = 1 TO M s = s + U(i, j) * B(i) NEXT i s = s / W(j) 'This is the divide by wj . END IF TMP(j) = s NEXT j FOR j = 1 TO N 'Matrix multiply by V to get answer. s = 0# FOR jj = 1 TO N s = s + V(j, jj) * TMP(jj) NEXT jj X(j) = s NEXT j RETURN 'Note that a typical use of svdcmp and svbksb superficially resembles the 'typical use of ludcmp and lubksb: In both cases, you decompose the left-hand 'matrix A just once, and then can use the decomposition either once or many times 'with different right-hand sides. The crucial difference is the "editing" of the singular 'values before using SVDCMP. 2000 'SUBROUTINE svdcmp(U,M,N,W,V) '-------------------------------------------------------------------------------------- ' Given a matrix A(1:m,1:n), this routine computes its singular value decomposition, ' A = U . W . Vt. The matrix U replaces A on output (1). The diagonal matrix of singular ' values W is output as a vector W(1:n). The matrix V (not the transpose Vt) is output ' as V(1:n,1:n). (1) here, matrix A is always called U. '-------------------------------------------------------------------------------------- DIM rv1(NMAX) g = 0# 'Householder reduction to bidiagonal form. scale = 0# anorm = 0# FOR i = 1 TO N l = i + 1 rv1(i) = scale * g g = 0# s = 0# scale = 0# IF i <= M THEN FOR k = i TO M scale = scale + ABS(U(k, i)) NEXT k IF scale <> 0# THEN FOR k = i TO M U(k, i) = U(k, i) / scale s = s + U(k, i) * U(k, i) NEXT k F = U(i, i) xa = SQR(s): xb = F: GOSUB 600 g = -Sign h = F * g - s U(i, i) = F - g FOR j = l TO N s = 0# FOR k = i TO M s = s + U(k, i) * U(k, j) NEXT k F = s / h FOR k = i TO M U(k, j) = U(k, j) + F * U(k, i) NEXT k NEXT j FOR k = i TO M U(k, i) = scale * U(k, i) NEXT k END IF END IF W(i) = scale * g g = 0# s = 0# scale = 0# IF i <= m.and.i <> N THEN FOR k = l TO N scale = scale + ABS(U(i, k)) NEXT k IF scale <> 0# THEN FOR k = l TO N U(i, k) = U(i, k) / scale s = s + U(i, k) * U(i, k) NEXT k F = U(i, l) xa = SQR(s): xb = F: GOSUB 600 g = -Sign h = F * g - s U(i, l) = F - g FOR k = l TO N rv1(k) = U(i, k) / h NEXT k FOR j = l TO M s = 0# FOR k = l TO N s = s + U(j, k) * U(i, k) NEXT k FOR k = l TO N U(j, k) = U(j, k) + s * rv1(k) NEXT k NEXT j FOR k = l TO N U(i, k) = scale * U(i, k) NEXT k END IF END IF 'anorm=Max(anorm,(abs(w(i))+abs(rv1(i)))) IF ABS(W(i)) >= ABS(rv1(i)) THEN anorm = ABS(W(i)) ELSE anorm = ABS(rv1(i)) END IF NEXT i FOR i = N TO 1 STEP -1'Accumulation of right-hand transformations. IF i < N THEN IF g <> 0# THEN FOR j = l TO N 'Double division to avoid possible underflow. V(j, i) = (U(i, j) / U(i, l)) / g NEXT j FOR j = l TO N s = 0# FOR k = l TO N s = s + U(i, k) * V(k, j) NEXT k FOR k = l TO N V(k, j) = V(k, j) + s * V(k, i) NEXT k NEXT j END IF FOR j = l TO N V(i, j) = 0# V(j, i) = 0# NEXT j END IF V(i, i) = 1# g = rv1(i) l = i NEXT i IF M < N THEN ii = M ELSE ii = N END IF FOR i = ii TO 1 STEP -1'Accumulation of left-hand transformations. l = i + 1 g = W(i) FOR j = l TO N U(i, j) = 0# NEXT j IF g <> 0# THEN g = 1# / g FOR j = l TO N s = 0# FOR k = l TO M s = s + U(k, i) * U(k, j) NEXT k F = (s / U(i, i)) * g FOR k = i TO M U(k, j) = U(k, j) + F * U(k, i) NEXT k NEXT j FOR j = i TO M U(j, i) = U(j, i) * g NEXT j ELSE FOR j = i TO M U(j, i) = 0# NEXT j END IF U(i, i) = U(i, i) + 1# NEXT i FOR k = N TO 1 STEP -1 'Diagonalization of the bidiagonal form: Loop over 'singular values, and over allowed iterations. FOR its = 1 TO 30 FOR l = k TO 1 STEP -1 'Test for splitting. nm = l - 1 'Note that rv1(1) is always zero. IF (ABS(rv1(l)) + anorm) = anorm THEN GOTO 2 IF (ABS(W(nm)) + anorm) = anorm THEN GOTO 1 NEXT l 1 c = 0# 'Cancellation of rv1(l), if l > 1. s = 1# FOR i = l TO k F = s * rv1(i) rv1(i) = c * rv1(i) IF (ABS(F) + anorm) = anorm THEN GOTO 2 g = W(i) xa = F: xb = g: GOSUB 500 h = Pythag W(i) = h h = 1# / h c = (g * h) s = -(F * h) FOR j = 1 TO M y = U(j, nm) z = U(j, i) U(j, nm) = (y * c) + (z * s) U(j, i) = -(y * s) + (z * c) NEXT j NEXT i 2 z = W(k) IF l = k THEN 'Convergence. IF z < 0# THEN 'Singular value is made nonnegative. W(k) = -z FOR j = 1 TO N V(j, k) = -V(j, k) NEXT j END IF GOTO 3 END IF IF (its.eq.30) THEN PRINT "No convergence in svdcmp." X = W(l) 'Shift from bottom 2-by-2 minor. nm = k - 1 y = W(nm) g = rv1(nm) h = rv1(k) F = ((y - z) * (y + z) + (g - h) * (g + h)) / (2! * h * y) xa = F: xb = 1#: GOSUB 500 g = Pythag xa = g: xb = F: GOSUB 600 F = ((X - z) * (X + z) + h * ((y / (F + Sign)) - h)) / X c = 1# 'Next QR transformation: s = 1# FOR j = l TO nm i = j + 1 g = rv1(i) y = W(i) h = s * g g = c * g xa = F: xb = h: GOSUB 500 z = Pythag 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 FOR jj = 1 TO N X = V(jj, j) z = V(jj, i) V(jj, j) = (X * c) + (z * s) V(jj, i) = -(X * s) + (z * c) NEXT jj xa = F: xb = h: GOSUB 500 z = Pythag W(j) = z 'Rotation can be arbitrary if z = 0. IF z <> 0# THEN z = 1# / z c = F * z s = h * z END IF F = (c * g) + (s * y) X = -(s * g) + (c * y) FOR jj = 1 TO M y = U(jj, j) z = U(jj, i) U(jj, j) = (y * c) + (z * s) U(jj, i) = -(y * s) + (z * c) NEXT jj NEXT j rv1(l) = 0# rv1(k) = F W(k) = X NEXT its 3 NEXT k RETURN 'end of file tsvbksb.bas