Doxygen Source Code Documentation
eis_svd.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | svd_ (integer *m, integer *n, integer *lda, doublereal *a, doublereal *w, logical *matu, integer *ldu, doublereal *u, logical *matv, integer *ldv, doublereal *v, integer *ierr, doublereal *rv1) |
Variables | |
doublereal | c_b47 = 1. |
Function Documentation
|
Definition at line 12 of file eis_svd.c. References a, abs, c_b47, d_sign(), i1, l, max, pythag_(), scale, and v. Referenced by svd_double().
00015 { 00016 /* System generated locals */ 00017 integer a_dim1, a_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, 00018 i__3; 00019 doublereal d__1, d__2, d__3, d__4; 00020 00021 /* Builtin functions */ 00022 double sqrt(doublereal), d_sign(doublereal *, doublereal *); 00023 00024 /* Local variables */ 00025 static doublereal c__, f, g, h__; 00026 static integer i__, j, k, l; 00027 static doublereal s, x, y, z__, scale; 00028 static integer i1, k1, l1, ii, kk, ll, mn; 00029 extern doublereal pythag_(doublereal *, doublereal *); 00030 static integer its; 00031 static doublereal tst1, tst2; 00032 00033 00034 00035 /* this subroutine is a translation of the algol procedure svd, */ 00036 /* num. math. 14, 403-420(1970) by golub and reinsch. */ 00037 /* handbook for auto. comp., vol ii-linear algebra, 134-151(1971). */ 00038 00039 /* this subroutine determines the singular value decomposition */ 00040 /* t */ 00041 /* a=usv of a real m by n rectangular matrix. householder */ 00042 /* bidiagonalization and a variant of the qr algorithm are used. */ 00043 00044 /* on input */ 00045 00046 /* nm must be set to the row dimension of two-dimensional */ 00047 /* array parameters as declared in the calling program */ 00048 /* dimension statement. note that nm must be at least */ 00049 /* as large as the maximum of m and n. */ 00050 00051 /* m is the number of rows of a (and u). */ 00052 00053 /* n is the number of columns of a, u, and v */ 00054 00055 /* a contains the rectangular input matrix to be decomposed. */ 00056 00057 /* matu should be set to .true. if the u matrix in the */ 00058 /* decomposition is desired, and to .false. otherwise. */ 00059 00060 /* matv should be set to .true. if the v matrix in the */ 00061 /* decomposition is desired, and to .false. otherwise. */ 00062 00063 /* lda, ldu, ldv: are the leading dimensions of matrices */ 00064 /* a, u, and v (respectively); must have */ 00065 /* lda >= m ; ldu >= m ; ldv >= n */ 00066 00067 /* on output */ 00068 00069 /* a is unaltered (unless overwritten by u or v). */ 00070 00071 /* w contains the n (non-negative) singular values of a (the */ 00072 /* diagonal elements of s). they are unordered. if an */ 00073 /* error exit is made, the singular values should be correct */ 00074 /* for indices ierr+1,ierr+2,...,n. */ 00075 00076 /* u contains the matrix u (orthogonal column vectors) of the */ 00077 /* decomposition if matu has been set to .true. otherwise */ 00078 /* u is used as a temporary array. u may coincide with a. */ 00079 /* if an error exit is made, the columns of u corresponding */ 00080 /* to indices of correct singular values should be correct. */ 00081 00082 /* v contains the matrix v (orthogonal) of the decomposition if */ 00083 /* matv has been set to .true. otherwise v is not referenced. */ 00084 /* v may also coincide with a if u is not needed. if an error */ 00085 /* exit is made, the columns of v corresponding to indices of */ 00086 /* correct singular values should be correct. */ 00087 00088 /* ierr is set to */ 00089 /* zero for normal return, */ 00090 /* k if the k-th singular value has not been */ 00091 /* determined after 30 iterations. */ 00092 00093 /* rv1 is a temporary storage array. */ 00094 00095 /* calls pythag for dsqrt(a*a + b*b) . */ 00096 00097 /* questions and comments should be directed to burton s. garbow, */ 00098 /* mathematics and computer science div, argonne national laboratory 00099 */ 00100 00101 /* this version dated august 1983. */ 00102 00103 /* ------------------------------------------------------------------ 00104 */ 00105 00106 /* Parameter adjustments */ 00107 --rv1; 00108 --w; 00109 a_dim1 = *lda; 00110 a_offset = a_dim1 + 1; 00111 a -= a_offset; 00112 u_dim1 = *ldu; 00113 u_offset = u_dim1 + 1; 00114 u -= u_offset; 00115 v_dim1 = *ldv; 00116 v_offset = v_dim1 + 1; 00117 if( v != (doublereal *)0 ) v -= v_offset; 00118 00119 /* Function Body */ 00120 *ierr = 0; 00121 00122 i__1 = *m; 00123 for (i__ = 1; i__ <= i__1; ++i__) { 00124 00125 i__2 = *n; 00126 for (j = 1; j <= i__2; ++j) { 00127 u[i__ + j * u_dim1] = a[i__ + j * a_dim1]; 00128 /* L100: */ 00129 } 00130 } 00131 /* .......... householder reduction to bidiagonal form .......... */ 00132 g = 0.; 00133 scale = 0.; 00134 x = 0.; 00135 00136 i__2 = *n; 00137 for (i__ = 1; i__ <= i__2; ++i__) { 00138 l = i__ + 1; 00139 rv1[i__] = scale * g; 00140 g = 0.; 00141 s = 0.; 00142 scale = 0.; 00143 if (i__ > *m) { 00144 goto L210; 00145 } 00146 00147 i__1 = *m; 00148 for (k = i__; k <= i__1; ++k) { 00149 /* L120: */ 00150 scale += (d__1 = u[k + i__ * u_dim1], abs(d__1)); 00151 } 00152 00153 if (scale == 0.) { 00154 goto L210; 00155 } 00156 00157 i__1 = *m; 00158 for (k = i__; k <= i__1; ++k) { 00159 u[k + i__ * u_dim1] /= scale; 00160 /* Computing 2nd power */ 00161 d__1 = u[k + i__ * u_dim1]; 00162 s += d__1 * d__1; 00163 /* L130: */ 00164 } 00165 00166 f = u[i__ + i__ * u_dim1]; 00167 d__1 = sqrt(s); 00168 g = -d_sign(&d__1, &f); 00169 h__ = f * g - s; 00170 u[i__ + i__ * u_dim1] = f - g; 00171 if (i__ == *n) { 00172 goto L190; 00173 } 00174 00175 i__1 = *n; 00176 for (j = l; j <= i__1; ++j) { 00177 s = 0.; 00178 00179 i__3 = *m; 00180 for (k = i__; k <= i__3; ++k) { 00181 /* L140: */ 00182 s += u[k + i__ * u_dim1] * u[k + j * u_dim1]; 00183 } 00184 00185 f = s / h__; 00186 00187 i__3 = *m; 00188 for (k = i__; k <= i__3; ++k) { 00189 u[k + j * u_dim1] += f * u[k + i__ * u_dim1]; 00190 /* L150: */ 00191 } 00192 } 00193 00194 L190: 00195 i__3 = *m; 00196 for (k = i__; k <= i__3; ++k) { 00197 /* L200: */ 00198 u[k + i__ * u_dim1] = scale * u[k + i__ * u_dim1]; 00199 } 00200 00201 L210: 00202 w[i__] = scale * g; 00203 g = 0.; 00204 s = 0.; 00205 scale = 0.; 00206 if (i__ > *m || i__ == *n) { 00207 goto L290; 00208 } 00209 00210 i__3 = *n; 00211 for (k = l; k <= i__3; ++k) { 00212 /* L220: */ 00213 scale += (d__1 = u[i__ + k * u_dim1], abs(d__1)); 00214 } 00215 00216 if (scale == 0.) { 00217 goto L290; 00218 } 00219 00220 i__3 = *n; 00221 for (k = l; k <= i__3; ++k) { 00222 u[i__ + k * u_dim1] /= scale; 00223 /* Computing 2nd power */ 00224 d__1 = u[i__ + k * u_dim1]; 00225 s += d__1 * d__1; 00226 /* L230: */ 00227 } 00228 00229 f = u[i__ + l * u_dim1]; 00230 d__1 = sqrt(s); 00231 g = -d_sign(&d__1, &f); 00232 h__ = f * g - s; 00233 u[i__ + l * u_dim1] = f - g; 00234 00235 i__3 = *n; 00236 for (k = l; k <= i__3; ++k) { 00237 /* L240: */ 00238 rv1[k] = u[i__ + k * u_dim1] / h__; 00239 } 00240 00241 if (i__ == *m) { 00242 goto L270; 00243 } 00244 00245 i__3 = *m; 00246 for (j = l; j <= i__3; ++j) { 00247 s = 0.; 00248 00249 i__1 = *n; 00250 for (k = l; k <= i__1; ++k) { 00251 /* L250: */ 00252 s += u[j + k * u_dim1] * u[i__ + k * u_dim1]; 00253 } 00254 00255 i__1 = *n; 00256 for (k = l; k <= i__1; ++k) { 00257 u[j + k * u_dim1] += s * rv1[k]; 00258 /* L260: */ 00259 } 00260 } 00261 00262 L270: 00263 i__1 = *n; 00264 for (k = l; k <= i__1; ++k) { 00265 /* L280: */ 00266 u[i__ + k * u_dim1] = scale * u[i__ + k * u_dim1]; 00267 } 00268 00269 L290: 00270 /* Computing MAX */ 00271 d__3 = x, d__4 = (d__1 = w[i__], abs(d__1)) + (d__2 = rv1[i__], abs( 00272 d__2)); 00273 x = max(d__3,d__4); 00274 /* L300: */ 00275 } 00276 /* .......... accumulation of right-hand transformations .......... */ 00277 if (! (*matv)) { 00278 goto L410; 00279 } 00280 /* .......... for i=n step -1 until 1 do -- .......... */ 00281 i__2 = *n; 00282 for (ii = 1; ii <= i__2; ++ii) { 00283 i__ = *n + 1 - ii; 00284 if (i__ == *n) { 00285 goto L390; 00286 } 00287 if (g == 0.) { 00288 goto L360; 00289 } 00290 00291 i__1 = *n; 00292 for (j = l; j <= i__1; ++j) { 00293 /* .......... double division avoids possible underflow ...... 00294 .... */ 00295 /* L320: */ 00296 v[j + i__ * v_dim1] = u[i__ + j * u_dim1] / u[i__ + l * u_dim1] / 00297 g; 00298 } 00299 00300 i__1 = *n; 00301 for (j = l; j <= i__1; ++j) { 00302 s = 0.; 00303 00304 i__3 = *n; 00305 for (k = l; k <= i__3; ++k) { 00306 /* L340: */ 00307 s += u[i__ + k * u_dim1] * v[k + j * v_dim1]; 00308 } 00309 00310 i__3 = *n; 00311 for (k = l; k <= i__3; ++k) { 00312 v[k + j * v_dim1] += s * v[k + i__ * v_dim1]; 00313 /* L350: */ 00314 } 00315 } 00316 00317 L360: 00318 i__3 = *n; 00319 for (j = l; j <= i__3; ++j) { 00320 v[i__ + j * v_dim1] = 0.; 00321 v[j + i__ * v_dim1] = 0.; 00322 /* L380: */ 00323 } 00324 00325 L390: 00326 v[i__ + i__ * v_dim1] = 1.; 00327 g = rv1[i__]; 00328 l = i__; 00329 /* L400: */ 00330 } 00331 /* .......... accumulation of left-hand transformations .......... */ 00332 L410: 00333 if (! (*matu)) { 00334 goto L510; 00335 } 00336 /* ..........for i=min(m,n) step -1 until 1 do -- .......... */ 00337 mn = *n; 00338 if (*m < *n) { 00339 mn = *m; 00340 } 00341 00342 i__2 = mn; 00343 for (ii = 1; ii <= i__2; ++ii) { 00344 i__ = mn + 1 - ii; 00345 l = i__ + 1; 00346 g = w[i__]; 00347 if (i__ == *n) { 00348 goto L430; 00349 } 00350 00351 i__3 = *n; 00352 for (j = l; j <= i__3; ++j) { 00353 /* L420: */ 00354 u[i__ + j * u_dim1] = 0.; 00355 } 00356 00357 L430: 00358 if (g == 0.) { 00359 goto L475; 00360 } 00361 if (i__ == mn) { 00362 goto L460; 00363 } 00364 00365 i__3 = *n; 00366 for (j = l; j <= i__3; ++j) { 00367 s = 0.; 00368 00369 i__1 = *m; 00370 for (k = l; k <= i__1; ++k) { 00371 /* L440: */ 00372 s += u[k + i__ * u_dim1] * u[k + j * u_dim1]; 00373 } 00374 /* .......... double division avoids possible underflow ...... 00375 .... */ 00376 f = s / u[i__ + i__ * u_dim1] / g; 00377 00378 i__1 = *m; 00379 for (k = i__; k <= i__1; ++k) { 00380 u[k + j * u_dim1] += f * u[k + i__ * u_dim1]; 00381 /* L450: */ 00382 } 00383 } 00384 00385 L460: 00386 i__1 = *m; 00387 for (j = i__; j <= i__1; ++j) { 00388 /* L470: */ 00389 u[j + i__ * u_dim1] /= g; 00390 } 00391 00392 goto L490; 00393 00394 L475: 00395 i__1 = *m; 00396 for (j = i__; j <= i__1; ++j) { 00397 /* L480: */ 00398 u[j + i__ * u_dim1] = 0.; 00399 } 00400 00401 L490: 00402 u[i__ + i__ * u_dim1] += 1.; 00403 /* L500: */ 00404 } 00405 /* .......... diagonalization of the bidiagonal form .......... */ 00406 L510: 00407 tst1 = x; 00408 /* .......... for k=n step -1 until 1 do -- .......... */ 00409 i__2 = *n; 00410 for (kk = 1; kk <= i__2; ++kk) { 00411 k1 = *n - kk; 00412 k = k1 + 1; 00413 its = 0; 00414 /* .......... test for splitting. */ 00415 /* for l=k step -1 until 1 do -- .......... */ 00416 L520: 00417 i__1 = k; 00418 for (ll = 1; ll <= i__1; ++ll) { 00419 l1 = k - ll; 00420 l = l1 + 1; 00421 tst2 = tst1 + (d__1 = rv1[l], abs(d__1)); 00422 if (tst2 == tst1) { 00423 goto L565; 00424 } 00425 /* .......... rv1(1) is always zero, so there is no exit */ 00426 /* through the bottom of the loop .......... */ 00427 tst2 = tst1 + (d__1 = w[l1], abs(d__1)); 00428 if (tst2 == tst1) { 00429 goto L540; 00430 } 00431 /* L530: */ 00432 } 00433 /* .......... cancellation of rv1(l) if l greater than 1 ......... 00434 . */ 00435 L540: 00436 c__ = 0.; 00437 s = 1.; 00438 00439 i__1 = k; 00440 for (i__ = l; i__ <= i__1; ++i__) { 00441 f = s * rv1[i__]; 00442 rv1[i__] = c__ * rv1[i__]; 00443 tst2 = tst1 + abs(f); 00444 if (tst2 == tst1) { 00445 goto L565; 00446 } 00447 g = w[i__]; 00448 h__ = pythag_(&f, &g); 00449 w[i__] = h__; 00450 c__ = g / h__; 00451 s = -f / h__; 00452 if (! (*matu)) { 00453 goto L560; 00454 } 00455 00456 i__3 = *m; 00457 for (j = 1; j <= i__3; ++j) { 00458 y = u[j + l1 * u_dim1]; 00459 z__ = u[j + i__ * u_dim1]; 00460 u[j + l1 * u_dim1] = y * c__ + z__ * s; 00461 u[j + i__ * u_dim1] = -y * s + z__ * c__; 00462 /* L550: */ 00463 } 00464 00465 L560: 00466 ; 00467 } 00468 /* .......... test for convergence .......... */ 00469 L565: 00470 z__ = w[k]; 00471 if (l == k) { 00472 goto L650; 00473 } 00474 /* .......... shift from bottom 2 by 2 minor .......... */ 00475 if (its == 30) { 00476 goto L1000; 00477 } 00478 ++its; 00479 x = w[l]; 00480 y = w[k1]; 00481 g = rv1[k1]; 00482 h__ = rv1[k]; 00483 f = ((g + z__) / h__ * ((g - z__) / y) + y / h__ - h__ / y) * .5; 00484 g = pythag_(&f, &c_b47); 00485 f = x - z__ / x * z__ + h__ / x * (y / (f + d_sign(&g, &f)) - h__); 00486 /* .......... next qr transformation .......... */ 00487 c__ = 1.; 00488 s = 1.; 00489 00490 i__1 = k1; 00491 for (i1 = l; i1 <= i__1; ++i1) { 00492 i__ = i1 + 1; 00493 g = rv1[i__]; 00494 y = w[i__]; 00495 h__ = s * g; 00496 g = c__ * g; 00497 z__ = pythag_(&f, &h__); 00498 rv1[i1] = z__; 00499 c__ = f / z__; 00500 s = h__ / z__; 00501 f = x * c__ + g * s; 00502 g = -x * s + g * c__; 00503 h__ = y * s; 00504 y *= c__; 00505 if (! (*matv)) { 00506 goto L575; 00507 } 00508 00509 i__3 = *n; 00510 for (j = 1; j <= i__3; ++j) { 00511 x = v[j + i1 * v_dim1]; 00512 z__ = v[j + i__ * v_dim1]; 00513 v[j + i1 * v_dim1] = x * c__ + z__ * s; 00514 v[j + i__ * v_dim1] = -x * s + z__ * c__; 00515 /* L570: */ 00516 } 00517 00518 L575: 00519 z__ = pythag_(&f, &h__); 00520 w[i1] = z__; 00521 /* .......... rotation can be arbitrary if z is zero ......... 00522 . */ 00523 if (z__ == 0.) { 00524 goto L580; 00525 } 00526 c__ = f / z__; 00527 s = h__ / z__; 00528 L580: 00529 f = c__ * g + s * y; 00530 x = -s * g + c__ * y; 00531 if (! (*matu)) { 00532 goto L600; 00533 } 00534 00535 i__3 = *m; 00536 for (j = 1; j <= i__3; ++j) { 00537 y = u[j + i1 * u_dim1]; 00538 z__ = u[j + i__ * u_dim1]; 00539 u[j + i1 * u_dim1] = y * c__ + z__ * s; 00540 u[j + i__ * u_dim1] = -y * s + z__ * c__; 00541 /* L590: */ 00542 } 00543 00544 L600: 00545 ; 00546 } 00547 00548 rv1[l] = 0.; 00549 rv1[k] = f; 00550 w[k] = x; 00551 goto L520; 00552 /* .......... convergence .......... */ 00553 L650: 00554 if (z__ >= 0.) { 00555 goto L700; 00556 } 00557 /* .......... w(k) is made non-negative .......... */ 00558 w[k] = -z__; 00559 if (! (*matv)) { 00560 goto L700; 00561 } 00562 00563 i__1 = *n; 00564 for (j = 1; j <= i__1; ++j) { 00565 /* L690: */ 00566 v[j + k * v_dim1] = -v[j + k * v_dim1]; 00567 } 00568 00569 L700: 00570 ; 00571 } 00572 00573 goto L1001; 00574 /* .......... set error -- no convergence to a */ 00575 /* singular value after 30 iterations .......... */ 00576 L1000: 00577 *ierr = k; 00578 L1001: 00579 return 0; 00580 } /* svd_ */ |
Variable Documentation
|
Definition at line 10 of file eis_svd.c. Referenced by svd_(). |