Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_svd.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008
00009
00010 static doublereal c_b47 = 1.;
00011
00012 int svd_(integer *m, integer *n, integer *lda, doublereal *a,
00013 doublereal *w, logical *matu, integer *ldu, doublereal *u, logical *
00014 matv, integer *ldv, doublereal *v, integer *ierr, doublereal *rv1)
00015 {
00016
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
00022 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00023
00024
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
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
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
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
00129 }
00130 }
00131
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
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
00161 d__1 = u[k + i__ * u_dim1];
00162 s += d__1 * d__1;
00163
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
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
00191 }
00192 }
00193
00194 L190:
00195 i__3 = *m;
00196 for (k = i__; k <= i__3; ++k) {
00197
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
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
00224 d__1 = u[i__ + k * u_dim1];
00225 s += d__1 * d__1;
00226
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
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
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
00259 }
00260 }
00261
00262 L270:
00263 i__1 = *n;
00264 for (k = l; k <= i__1; ++k) {
00265
00266 u[i__ + k * u_dim1] = scale * u[i__ + k * u_dim1];
00267 }
00268
00269 L290:
00270
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
00275 }
00276
00277 if (! (*matv)) {
00278 goto L410;
00279 }
00280
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
00294
00295
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
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
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
00323 }
00324
00325 L390:
00326 v[i__ + i__ * v_dim1] = 1.;
00327 g = rv1[i__];
00328 l = i__;
00329
00330 }
00331
00332 L410:
00333 if (! (*matu)) {
00334 goto L510;
00335 }
00336
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
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
00372 s += u[k + i__ * u_dim1] * u[k + j * u_dim1];
00373 }
00374
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
00382 }
00383 }
00384
00385 L460:
00386 i__1 = *m;
00387 for (j = i__; j <= i__1; ++j) {
00388
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
00398 u[j + i__ * u_dim1] = 0.;
00399 }
00400
00401 L490:
00402 u[i__ + i__ * u_dim1] += 1.;
00403
00404 }
00405
00406 L510:
00407 tst1 = x;
00408
00409 i__2 = *n;
00410 for (kk = 1; kk <= i__2; ++kk) {
00411 k1 = *n - kk;
00412 k = k1 + 1;
00413 its = 0;
00414
00415
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
00426
00427 tst2 = tst1 + (d__1 = w[l1], abs(d__1));
00428 if (tst2 == tst1) {
00429 goto L540;
00430 }
00431
00432 }
00433
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
00463 }
00464
00465 L560:
00466 ;
00467 }
00468
00469 L565:
00470 z__ = w[k];
00471 if (l == k) {
00472 goto L650;
00473 }
00474
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
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
00516 }
00517
00518 L575:
00519 z__ = pythag_(&f, &h__);
00520 w[i1] = z__;
00521
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
00542 }
00543
00544 L600:
00545 ;
00546 }
00547
00548 rv1[l] = 0.;
00549 rv1[k] = f;
00550 w[k] = x;
00551 goto L520;
00552
00553 L650:
00554 if (z__ >= 0.) {
00555 goto L700;
00556 }
00557
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
00566 v[j + k * v_dim1] = -v[j + k * v_dim1];
00567 }
00568
00569 L700:
00570 ;
00571 }
00572
00573 goto L1001;
00574
00575
00576 L1000:
00577 *ierr = k;
00578 L1001:
00579 return 0;
00580 }
00581