#version 3.8; global_settings {assumed_gamma 1.0} // Adapted to POV-Ray SDL from EISPACK Fortran90 code // Bill "Bald Eagle" Walker - April 2021 // https://www.netlib.org/eispack/svd.f #declare Sign = function (_A, _B) {select (_B, -abs(_A), abs(_A), abs(_A))} #declare pythag = function (_A, _B) {sqrt(pow(_A,2) + pow(_B,2))} /* subroutine svd (nm,m,n,a,w,matu,u,matv,v,ierr,rv1) integer i,j,k,l,m,n,ii,i1,kk,k1,ll,l1,mn,nm,its,ierr double precision a(nm,n),w(n),u(nm,n),V(nm,n),rv1(n) double precision c,f,g,h,s,x,y,z,tst1,tst2,scale,pythag logical matu,matv nm must be set to the row dimension of two-dimensional array parameters as declared in the calling program dimension statement. note that nm must be at least as large as the maximum of m and n. m is the number of rows of a (and u). n is the number of columns of a (and u) and the order of v. a contains the rectangular input matrix to be decomposed. matu should be set to .true. if the u matrix in the decomposition is desired, and to .false. otherwise. matv should be set to .true. if the v matrix in the decomposition is desired, and to .false. otherwise. rv1 is a temporary storage array. */ #macro Display2DArray (_Array, _Mname, _P) // display contents of vector array to _P decimal places #debug concat ("Matrix ", _Mname, " = \n") #debug "------------------------ \n" #local _rows = dimension_size (_Array, 1)-1; #local _columns = dimension_size (_Array, 2)-1; #for (_X, 0, _columns) #debug concat ( #for (_Y, 0, _rows) "[", str(_Array [_Y][_X], -(_P+2), _P), "]", #end "\n") #end #debug "------------------------ \n" #end #declare Test = array [4][3]{ {1, 2, 4}, {6, -1, 3}, {4, 0, -2}, {5, 2, 3} } #declare M = dimension_size (Test, 1); // A rows #declare N = dimension_size (Test, 2); // A columns #declare NM = max (M, N); Display2DArray (Test, "Test", 6) //#error "end" #macro Make2DArray (_N, _M, _val) #declare Mat = array [_M][_N]; #for (_col, 0, _N-1) #for (_row, 0, _M-1) #declare Mat[_row][_col] = _val; #end #end Mat #end #macro SVD (nm, m, n, A, matu, matv) #local rv1 = array; #local ierr = 0; // copy Matrix A into Matrix U #local Ai = dimension_size (A, 1)-1; #local Aj = dimension_size (A, 2)-1; #local U = Make2DArray (Ai+2, Ai+2, 0); #local V = Make2DArray (Ai+2, Aj+2, 0); /* #for (i, 0, Ai) #for (j, 0, Aj) #local U [i][j] = 0; #local V [i][j] = 0; #end // end for (j, 0, Aj) #end // end for (i, 0, Ai) */ #local w = array; //#for (i, 1, m) // do 100 // #for (j, 1, n) // do 100 // U [i, j] = A [i,j]; // #end //#end // 100 // .......... householder reduction to bidiagonal form .......... #local g = 0; #local _scale = 0; #local _x = 0; #for (i, 1, n) // do 300 #local l = i + 1; #local rv1 [i] = _scale * g; #local g = 0; #local s = 0; #local _scale = 0; #if (i > m) // go to 210 #else #for (k, i, m) // do 120 #local _scale = _scale + abs (U [k][i]); #if (_scale = 0) // go to 210 #else #for (k, i, m) // do 130 #local U [k][i] = U [k][i] / _scale; #local s = s + pow(U [k][i], 2); #end // end for (k = i, m) // 130 #local f = U [i][i]; #local g = -Sign (sqrt (s), f); #local h = f * g - s; #local U [i][i] = f - g; #if (i = n) // goto 190 #else #for (j, l, n) // do 150 #local s = 0; #for (k, i, m) // do 140 #local s = s + U [k][i] * U [k][j]; #end // end for // 140 #local f = s / h; #for (k, i, m) // do 150 #local U [k][j] = U [k][j] + f * U [k][i]; #end // end for #end // end for (j = l, n) // 150 #for (k, i, m) // 190 do 200 #local U [k][i] = _scale * U [k][i]; #end // end for (k = i, m) // 200 #end // end if (i = n) #end // end if (_scale = 0) // 210 #end // end for (k = i, m) #end // end if (i > m) // 210 #local w [i] = _scale * g; // 210 #local g = 0; #local s = 0; #local _scale = 0; #if (i > m | i = n) // go to 290 #else #for (k, l, n) // do 220 #local _scale = _scale + abs (U [i][k]); #end // 220 #if (_scale = 0) // go to 290 #else #for (k, l, n) // do 230 #local U [i][k] = U [i][k] / _scale; #local s = s + pow(U [i][k], 2); #end // 230 #local f = U [i][l]; #local g = -Sign (sqrt (s), f); #local h = f * g - s; #local U [i][l] = f - g; #for (k, l, n) // do 240 #local rv1(k) = U [i][k] / h; #end // 240 #if (i = m) // go to 270 #else #for (j, l, m) // do 260 #local s = 0; #for (k, l, n) // do 250 #local s = s + U [j][k] * U [i][k]; #end // 250 #for (k, l, n) // do 260 #local U [j][k] = U [j][k] + s * rv1 [k]; #end // 260 #end // 260 #end // 270 #for (k, l, n) // 270 do 280 #local U [i][k] = _scale * U [i][k]; #end // 280 #end // 290 #end //290 #local _x = max(_x, abs (w [i]) + abs (rv1 [i])); // 290 dmax1 #end // 300 // .......... accumulation of right-hand transformations .......... #if (!matv) // go to 410 #else // .......... for i=n step -1 until 1 do -- .......... #for (ii, 1, n) // do 400 #local i = n + 1 - ii; #if (i = n) // go to 390 #else #if (g = 0) // go to 360 #else #for (j, l, n) // do 320 // .......... double division avoids possible underflow .......... #local V [j][i] = (U [i][j] / U [i][l]) / g; #end // 320 #for (j, l, n) // do 350 #local s = 0; #for (k, l, n) // do 340 #local s = s + U [i][k] * V [k][j]; #end // 340 #for (k, l, n) // do 350 #local V [k][j] = V [k][j] + s * V [k][i]; #end // 350 #end // 350 #end // 360 #for (j, l, n) // 360 do 380 #local V [i][j] = 0; #local V [j][i] = 0; #end // 380 #end // 390 #local V [i][i] = 1; // 390 #local g = rv1 [i]; #local l = i; #end // 400 // .......... accumulation of left-hand transformations .......... #end // 410 #if (!matu) // 410 // go to 510 #else // ..........for i=min(m,n) step -1 until 1 do -- .......... #local mn = n; #if (m < n) #local mn = m; #end #for (ii, 1, mn) // do 500 #local i = mn + 1 - ii; #local l = i + 1; #local g = w [i]; #if (i = n) // go to 430 #end #else #for (j, l, n) // do 420 #local U [i][j] = 0; #end // 420 #end // 430 #if (g = 0) // 430 // go to 475 #else #if (i = mn) //go to 460 #end #else #for (j, l, n) // do 450 #local s = 0; #for (k, l, m) // do 440 #local s = s + U [k][i] * U [k][j]; #end // end for (k = l, m) // 440 // .......... double division avoids possible underflow .......... #local f = (s / U [i][i]) / g; #for (k, i, m) // do 450 #local U [k][j] = U [k][j] + f * U [k][i]; #end // end for (k = i, m) // 450 #end // end for(j = l, n) // 450 #end // end if (i = mn) // 460 #for (j, i, m) // 460 do 470 #local U [j][i] = U [j][i] / g; #end // 470 // go to 490 [BW: add secondary #if block to skip 475] #end // end if (g = 0) // 475 #if (i = mn) // BW #for (j, i, m) // 475 do 480 #local U [j][i] = 0; #end // 480 #end // BW #local U [i][i] = U [i][i] + 1; // 490 #end // 500 // .......... diagonalization of the bidiagonal form .......... #end // 510 #local tst1 = _x; // 510 // .......... for k=n step -1 until 1 do -- .......... #for (kk, 1, n) // do 700 #local k1 = n - kk; #local k = k1 + 1; #local its = 0; // .......... test for splitting. // for l=k step -1 until 1 do -- .......... #local Goto520 = 1; #while (Goto520 = 1) // BW #local Goto520 = 0; #for (ll, 1, k) // 520 do 530 #local l1 = k - ll; #local l = l1 + 1; #local tst2 = tst1 + abs (rv1 [l]); #if (tst2 = tst1) /* go to 565 */ #local BreakEarly = 1; #break #end // HMMMMMMMMMMMMMMMMMMMM // .......... rv1[1] is always zero, so there is no exit // through the bottom of the loop .......... #local tst2 = tst1 + abs (w [l1]); #if (tst2 = tst1) /* go to 540 (exit loop) */ #break #end //----------------------- #end // 530 // .......... cancellation of rv1(l) if l greater than 1 .......... #if (BreakEarly) // BW #else // BW #local c = 0; // 540 #local s = 1; #end // BW #for (i, l, k) // do 560 #if (BreakEarly) #break #end // (skip this loop and goto 565) BW #local f = s * rv1 [i]; #local rv1 [i] = c * rv1 [i]; #local tst2 = tst1 + abs (f); #if (tst2 = tst1) /* go to 565 (exit loop) */ #break #end // ----------------------- #local g = w [i]; #local h = pythag (f, g); #local w [i] = h; #local c = g / h; #local s = -f / h; #if (!matu) /* go to 560 (exit loop) */ #break #end #for (j, 1, m) // do 550 #local _y = U [j][l1]; #local _z = U [j][i]; #local U [j][l1] = _y * c + _z * s; #local U [j][i] = -_y * s + _z * c; #end // 550 #end // 560 // .......... test for convergence .......... #local _z = w[k]; // 565 #if (l = k) // -------------------- // go to 650 #else // .......... shift from bottom 2 by 2 minor .......... #if (its = 30) /* go to 1000 */ #break #end // -------------------- #local its = its + 1; #local _x = w [l]; #local _y = w [k1]; #local g = rv1 [k1]; #local h = rv1 [k]; #local f = 0.5 * (((g + _z) / h) * ((g - _z) / _y) + _y / h - h / _y); #local g = pythag (f, 1); #local f = _x - (_z / _x) * _z + (h / _x) * (_y / (f + Sign (g,f)) - h); // .......... next qr transformation .......... #local c = 1 #local s = 1 #for(i1, l, k1) // do 600 #local i = i1 + 1; #local g = rv1 [i]; #local _y = w [i]; #local h = s * g; #local g = c * g; #local _z = pythag (f, h); #local rv1 [i1] = _z; #local c = f / _z; #local s = h / _z; #local f = _x * c + g * s; #local g = -_x * s + g * c; #local h = _y * s; #local _y = _y * c; #if (!matv) // -------------------- // go to 575 #end #else #for (j, 1, n) // do 570 #local _x = V [j][i1]; #local _z = V [j][i]; #local V [j][i1] = _x * c + _z * s; #local V [j][i] = -_x * s + _z * c; #end // 570 #end // 575 #local _z = pythag (f, h); // 575 #local w [i1] = _z; // .......... rotation can be arbitrary if z is zero .......... #if (_z = 0) // go to 580 #else #local c = f / _z; #local s = h / _z; #end // 580 #local f = c * g + s * _y; // 580 #local _x = -s * g + c * _y; #if (!matu) // -------------------- // go to 600 #else #for (j, 1, m) // do 590 #local _y = U [j][i1] #local _z = U [j][i] #local U [j][i1] = _y * c + _z * s #local U [j][i] = -_y * s + _z * c #end // 590 #end // 600 #end // 600 #local rv1 [l] = 0; #local rv1 [k] = f; #local w [k] = _x; #local Goto520 = 1; //go to 520 // **** UH-Oh! **** // .......... convergence .......... #end // 650 #end // end while Goto520 BW #if (_z >= 0) // 650 -------------------- // go to 700 #else #if (its = 30) /* go to 1000 */ #break #end // BW // .......... w [k) is made non-negative .......... #local w [k] = -_z #if (!matv) // -------------------- // go to 700 #else #for (j, 1, n) // do 690 #local V [j][k] = -V [j][k]; #end #end // 700 #end // 700 #end // 700 //go to 1001 // -------------------- #if (its = 30) /* go to 1000 */ // BW // .......... set error -- no convergence to a // singular value after 30 iterations .......... #local ierr = k; // 1000 #else // BW // 1001 return #end // BW #declare UU = U; #declare WW = w; #declare VV = V; #end SVD (NM, M, N, Test, yes, yes) Display2DArray (UU, "U", 6) Display2DArray (VV, "V", 6) #error "end" /* subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1) c integer i,j,k,l,m,n,ii,i1,kk,k1,ll,l1,mn,nm,its,ierr double precision a(nm,n),w(n),u(nm,n),V(nm,n),rv1(n) double precision c,f,g,h,s,x,y,z,tst1,tst2,scale,pythag logical matu,matv c c this subroutine is a translation of the algol procedure svd, c num. math. 14, 403-420(1970) by golub and reinsch. c handbook for auto. comp., vol ii-linear algebra, 134-151(1971). c c this subroutine determines the singular value decomposition c t c a=usv of a real m by n rectangular matrix. householder c bidiagonalization and a variant of the qr algorithm are used. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. note that nm must be at least c as large as the maximum of m and n. c c m is the number of rows of a (and u). c c n is the number of columns of a (and u) and the order of v. c c a contains the rectangular input matrix to be decomposed. c c matu should be set to .true. if the u matrix in the c decomposition is desired, and to .false. otherwise. c c matv should be set to .true. if the v matrix in the c decomposition is desired, and to .false. otherwise. c c on output c c a is unaltered (unless overwritten by u or v). c c w contains the n (non-negative) singular values of a (the c diagonal elements of s). they are unordered. if an c error exit is made, the singular values should be correct c for indices ierr+1,ierr+2,...,n. c c u contains the matrix u (orthogonal column vectors) of the c decomposition if matu has been set to .true. otherwise c u is used as a temporary array. u may coincide with a. c if an error exit is made, the columns of u corresponding c to indices of correct singular values should be correct. c c v contains the matrix v (orthogonal) of the decomposition if c matv has been set to .true. otherwise v is not referenced. c v may also coincide with a if u is not needed. if an error c exit is made, the columns of v corresponding to indices of c correct singular values should be correct. c c ierr is set to c zero for normal return, c k if the k-th singular value has not been c determined after 30 iterations. c c rv1 is a temporary storage array. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 c do 100 i = 1, m c do 100 j = 1, n u(i,j) = a(i,j) 100 continue c .......... householder reduction to bidiagonal form .......... g = 0.0d0 scale = 0.0d0 x = 0.0d0 c do 300 i = 1, n l = i + 1 rv1(i) = scale * g g = 0.0d0 s = 0.0d0 scale = 0.0d0 if (i .gt. m) go to 210 c do 120 k = i, m 120 scale = scale + dabs(u(k,i)) c if (scale .eq. 0.0d0) go to 210 c do 130 k = i, m u(k,i) = u(k,i) / scale s = s + u(k,i)**2 130 continue c f = u(i,i) g = -dsign(dsqrt(s),f) h = f * g - s u(i,i) = f - g if (i .eq. n) go to 190 c do 150 j = l, n s = 0.0d0 c do 140 k = i, m 140 s = s + u(k,i) * u(k,j) c f = s / h c do 150 k = i, m u(k,j) = u(k,j) + f * u(k,i) 150 continue c 190 do 200 k = i, m 200 u(k,i) = scale * u(k,i) c 210 w(i) = scale * g g = 0.0d0 s = 0.0d0 scale = 0.0d0 if (i .gt. m .or. i .eq. n) go to 290 c do 220 k = l, n 220 scale = scale + dabs(u(i,k)) c if (scale .eq. 0.0d0) go to 290 c do 230 k = l, n u(i,k) = u(i,k) / scale s = s + u(i,k)**2 230 continue c f = u(i,l) g = -dsign(dsqrt(s),f) h = f * g - s u(i,l) = f - g c do 240 k = l, n 240 rv1(k) = u(i,k) / h c if (i .eq. m) go to 270 c do 260 j = l, m s = 0.0d0 c do 250 k = l, n 250 s = s + u(j,k) * u(i,k) c do 260 k = l, n u(j,k) = u(j,k) + s * rv1(k) 260 continue c 270 do 280 k = l, n 280 u(i,k) = scale * u(i,k) c 290 x = dmax1(x,dabs(w(i))+dabs(rv1(i))) 300 continue c .......... accumulation of right-hand transformations .......... if (.not. matv) go to 410 c .......... for i=n step -1 until 1 do -- .......... do 400 ii = 1, n i = n + 1 - ii if (i .eq. n) go to 390 if (g .eq. 0.0d0) go to 360 c do 320 j = l, n c .......... double division avoids possible underflow .......... 320 v(j,i) = (u(i,j) / u(i,l)) / g c do 350 j = l, n s = 0.0d0 c do 340 k = l, n 340 s = s + u(i,k) * v(k,j) c do 350 k = l, n v(k,j) = v(k,j) + s * v(k,i) 350 continue c 360 do 380 j = l, n v(i,j) = 0.0d0 v(j,i) = 0.0d0 380 continue c 390 v(i,i) = 1.0d0 g = rv1(i) l = i 400 continue c .......... accumulation of left-hand transformations .......... 410 if (.not. matu) go to 510 c ..........for i=min(m,n) step -1 until 1 do -- .......... mn = n if (m .lt. n) mn = m c do 500 ii = 1, mn i = mn + 1 - ii l = i + 1 g = w(i) if (i .eq. n) go to 430 c do 420 j = l, n 420 u(i,j) = 0.0d0 c 430 if (g .eq. 0.0d0) go to 475 if (i .eq. mn) go to 460 c do 450 j = l, n s = 0.0d0 c do 440 k = l, m 440 s = s + u(k,i) * u(k,j) c .......... double division avoids possible underflow .......... f = (s / u(i,i)) / g c do 450 k = i, m u(k,j) = u(k,j) + f * u(k,i) 450 continue c 460 do 470 j = i, m 470 u(j,i) = u(j,i) / g c go to 490 c 475 do 480 j = i, m 480 u(j,i) = 0.0d0 c 490 u(i,i) = u(i,i) + 1.0d0 500 continue c .......... diagonalization of the bidiagonal form .......... 510 tst1 = x c .......... for k=n step -1 until 1 do -- .......... do 700 kk = 1, n k1 = n - kk k = k1 + 1 its = 0 c .......... test for splitting. c for l=k step -1 until 1 do -- .......... 520 do 530 ll = 1, k l1 = k - ll l = l1 + 1 tst2 = tst1 + dabs(rv1(l)) if (tst2 .eq. tst1) go to 565 c .......... rv1(1) is always zero, so there is no exit c through the bottom of the loop .......... tst2 = tst1 + dabs(w(l1)) if (tst2 .eq. tst1) go to 540 530 continue c .......... cancellation of rv1(l) if l greater than 1 .......... 540 c = 0.0d0 s = 1.0d0 c do 560 i = l, k f = s * rv1(i) rv1(i) = c * rv1(i) tst2 = tst1 + dabs(f) if (tst2 .eq. tst1) go to 565 g = w(i) h = pythag(f,g) w(i) = h c = g / h s = -f / h if (.not. matu) go to 560 c do 550 j = 1, m y = u(j,l1) z = u(j,i) u(j,l1) = y * c + z * s u(j,i) = -y * s + z * c 550 continue c 560 continue c .......... test for convergence .......... 565 z = w(k) if (l .eq. k) go to 650 c .......... shift from bottom 2 by 2 minor .......... if (its .eq. 30) go to 1000 its = its + 1 x = w(l) y = w(k1) g = rv1(k1) h = rv1(k) f = 0.5d0 * (((g + z) / h) * ((g - z) / y) + y / h - h / y) g = pythag(f,1.0d0) f = x - (z / x) * z + (h / x) * (y / (f + dsign(g,f)) - h) c .......... next qr transformation .......... c = 1.0d0 s = 1.0d0 c do 600 i1 = l, k1 i = i1 + 1 g = rv1(i) y = w(i) h = s * g g = c * g z = pythag(f,h) rv1(i1) = z c = f / z s = h / z f = x * c + g * s g = -x * s + g * c h = y * s y = y * c if (.not. matv) go to 575 c do 570 j = 1, n x = v(j,i1) z = v(j,i) v(j,i1) = x * c + z * s v(j,i) = -x * s + z * c 570 continue c 575 z = pythag(f,h) w(i1) = z c .......... rotation can be arbitrary if z is zero .......... if (z .eq. 0.0d0) go to 580 c = f / z s = h / z 580 f = c * g + s * y x = -s * g + c * y if (.not. matu) go to 600 c do 590 j = 1, m y = u(j,i1) z = u(j,i) u(j,i1) = y * c + z * s u(j,i) = -y * s + z * c 590 continue c 600 continue c rv1(l) = 0.0d0 rv1(k) = f w(k) = x go to 520 c .......... convergence .......... 650 if (z .ge. 0.0d0) go to 700 c .......... w(k) is made non-negative .......... w(k) = -z if (.not. matv) go to 700 c do 690 j = 1, n 690 v(j,k) = -v(j,k) c 700 continue c go to 1001 c .......... set error -- no convergence to a c singular value after 30 iterations .......... 1000 ierr = k 1001 return end */ #debug "\n"