📄 csvdc.c
字号:
, dabs(r__2)) == (float)0.) {
goto L120;
}
/* apply the transformation. */
/*< do 90 i = lp1, n >*/
i__2 = *n;
for (i__ = lp1; i__ <= i__2; ++i__) {
/*< work(i) = (0.0e0,0.0e0) >*/
i__3 = i__;
work[i__3].r = (float)0., work[i__3].i = (float)0.;
/*< 90 continue >*/
/* L90: */
}
/*< do 100 j = lp1, p >*/
i__2 = *p;
for (j = lp1; j <= i__2; ++j) {
/*< call caxpy(n-l,e(j),x(lp1,j),1,work(lp1),1) >*/
i__3 = *n - l;
caxpy_(&i__3, &e[j], &x[lp1 + j * x_dim1], &c__1, &work[lp1], &
c__1);
/*< 100 continue >*/
/* L100: */
}
/*< do 110 j = lp1, p >*/
i__2 = *p;
for (j = lp1; j <= i__2; ++j) {
/*< >*/
i__3 = *n - l;
i__4 = j;
q__3.r = -e[i__4].r, q__3.i = -e[i__4].i;
c_div(&q__2, &q__3, &e[lp1]);
r_cnjg(&q__1, &q__2);
caxpy_(&i__3, &q__1, &work[lp1], &c__1, &x[lp1 + j * x_dim1], &
c__1);
/*< 110 continue >*/
/* L110: */
}
/*< 120 continue >*/
L120:
/*< if (.not.wantv) go to 140 >*/
if (! wantv) {
goto L140;
}
/* place the transformation in v for subsequent */
/* back multiplication. */
/*< do 130 i = lp1, p >*/
i__2 = *p;
for (i__ = lp1; i__ <= i__2; ++i__) {
/*< v(i,l) = e(i) >*/
i__3 = i__ + l * v_dim1;
i__4 = i__;
v[i__3].r = e[i__4].r, v[i__3].i = e[i__4].i;
/*< 130 continue >*/
/* L130: */
}
/*< 140 continue >*/
L140:
/*< 150 continue >*/
L150:
/*< 160 continue >*/
/* L160: */
;
}
/*< 170 continue >*/
L170:
/* set up the final bidiagonal matrix or order m. */
/*< m = min0(p,n+1) >*/
/* Computing MIN */
i__1 = *p, i__2 = *n + 1;
m = min(i__1,i__2);
/*< nctp1 = nct + 1 >*/
nctp1 = nct + 1;
/*< nrtp1 = nrt + 1 >*/
nrtp1 = nrt + 1;
/*< if (nct .lt. p) s(nctp1) = x(nctp1,nctp1) >*/
if (nct < *p) {
i__1 = nctp1;
i__2 = nctp1 + nctp1 * x_dim1;
s[i__1].r = x[i__2].r, s[i__1].i = x[i__2].i;
}
/*< if (n .lt. m) s(m) = (0.0e0,0.0e0) >*/
if (*n < m) {
i__1 = m;
s[i__1].r = (float)0., s[i__1].i = (float)0.;
}
/*< if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m) >*/
if (nrtp1 < m) {
i__1 = nrtp1;
i__2 = nrtp1 + m * x_dim1;
e[i__1].r = x[i__2].r, e[i__1].i = x[i__2].i;
}
/*< e(m) = (0.0e0,0.0e0) >*/
i__1 = m;
e[i__1].r = (float)0., e[i__1].i = (float)0.;
/* if required, generate u. */
/*< if (.not.wantu) go to 300 >*/
if (! wantu) {
goto L300;
}
/*< if (ncu .lt. nctp1) go to 200 >*/
if (ncu < nctp1) {
goto L200;
}
/*< do 190 j = nctp1, ncu >*/
i__1 = ncu;
for (j = nctp1; j <= i__1; ++j) {
/*< do 180 i = 1, n >*/
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
/*< u(i,j) = (0.0e0,0.0e0) >*/
i__3 = i__ + j * u_dim1;
u[i__3].r = (float)0., u[i__3].i = (float)0.;
/*< 180 continue >*/
/* L180: */
}
/*< u(j,j) = (1.0e0,0.0e0) >*/
i__2 = j + j * u_dim1;
u[i__2].r = (float)1., u[i__2].i = (float)0.;
/*< 190 continue >*/
/* L190: */
}
/*< 200 continue >*/
L200:
/*< if (nct .lt. 1) go to 290 >*/
if (nct < 1) {
goto L290;
}
/*< do 280 ll = 1, nct >*/
i__1 = nct;
for (ll = 1; ll <= i__1; ++ll) {
/*< l = nct - ll + 1 >*/
l = nct - ll + 1;
/*< if (cabs1(s(l)) .eq. 0.0e0) go to 250 >*/
i__2 = l;
if ((r__1 = s[i__2].r, dabs(r__1)) + (r__2 = r_imag(&s[l]), dabs(r__2)
) == (float)0.) {
goto L250;
}
/*< lp1 = l + 1 >*/
lp1 = l + 1;
/*< if (ncu .lt. lp1) go to 220 >*/
if (ncu < lp1) {
goto L220;
}
/*< do 210 j = lp1, ncu >*/
i__2 = ncu;
for (j = lp1; j <= i__2; ++j) {
/*< t = -cdotc(n-l+1,u(l,l),1,u(l,j),1)/u(l,l) >*/
i__3 = *n - l + 1;
cdotc_(&q__3, &i__3, &u[l + l * u_dim1], &c__1, &u[l + j * u_dim1]
, &c__1);
q__2.r = -q__3.r, q__2.i = -q__3.i;
c_div(&q__1, &q__2, &u[l + l * u_dim1]);
t.r = q__1.r, t.i = q__1.i;
/*< call caxpy(n-l+1,t,u(l,l),1,u(l,j),1) >*/
i__3 = *n - l + 1;
caxpy_(&i__3, &t, &u[l + l * u_dim1], &c__1, &u[l + j * u_dim1], &
c__1);
/*< 210 continue >*/
/* L210: */
}
/*< 220 continue >*/
L220:
/*< call cscal(n-l+1,(-1.0e0,0.0e0),u(l,l),1) >*/
i__2 = *n - l + 1;
cscal_(&i__2, &c_b53, &u[l + l * u_dim1], &c__1);
/*< u(l,l) = (1.0e0,0.0e0) + u(l,l) >*/
i__2 = l + l * u_dim1;
i__3 = l + l * u_dim1;
q__1.r = u[i__3].r + (float)1., q__1.i = u[i__3].i + (float)0.;
u[i__2].r = q__1.r, u[i__2].i = q__1.i;
/*< lm1 = l - 1 >*/
lm1 = l - 1;
/*< if (lm1 .lt. 1) go to 240 >*/
if (lm1 < 1) {
goto L240;
}
/*< do 230 i = 1, lm1 >*/
i__2 = lm1;
for (i__ = 1; i__ <= i__2; ++i__) {
/*< u(i,l) = (0.0e0,0.0e0) >*/
i__3 = i__ + l * u_dim1;
u[i__3].r = (float)0., u[i__3].i = (float)0.;
/*< 230 continue >*/
/* L230: */
}
/*< 240 continue >*/
L240:
/*< go to 270 >*/
goto L270;
/*< 250 continue >*/
L250:
/*< do 260 i = 1, n >*/
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
/*< u(i,l) = (0.0e0,0.0e0) >*/
i__3 = i__ + l * u_dim1;
u[i__3].r = (float)0., u[i__3].i = (float)0.;
/*< 260 continue >*/
/* L260: */
}
/*< u(l,l) = (1.0e0,0.0e0) >*/
i__2 = l + l * u_dim1;
u[i__2].r = (float)1., u[i__2].i = (float)0.;
/*< 270 continue >*/
L270:
/*< 280 continue >*/
/* L280: */
;
}
/*< 290 continue >*/
L290:
/*< 300 continue >*/
L300:
/* if it is required, generate v. */
/*< if (.not.wantv) go to 350 >*/
if (! wantv) {
goto L350;
}
/*< do 340 ll = 1, p >*/
i__1 = *p;
for (ll = 1; ll <= i__1; ++ll) {
/*< l = p - ll + 1 >*/
l = *p - ll + 1;
/*< lp1 = l + 1 >*/
lp1 = l + 1;
/*< if (l .gt. nrt) go to 320 >*/
if (l > nrt) {
goto L320;
}
/*< if (cabs1(e(l)) .eq. 0.0e0) go to 320 >*/
i__2 = l;
if ((r__1 = e[i__2].r, dabs(r__1)) + (r__2 = r_imag(&e[l]), dabs(r__2)
) == (float)0.) {
goto L320;
}
/*< do 310 j = lp1, p >*/
i__2 = *p;
for (j = lp1; j <= i__2; ++j) {
/*< t = -cdotc(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l) >*/
i__3 = *p - l;
cdotc_(&q__3, &i__3, &v[lp1 + l * v_dim1], &c__1, &v[lp1 + j *
v_dim1], &c__1);
q__2.r = -q__3.r, q__2.i = -q__3.i;
c_div(&q__1, &q__2, &v[lp1 + l * v_dim1]);
t.r = q__1.r, t.i = q__1.i;
/*< call caxpy(p-l,t,v(lp1,l),1,v(lp1,j),1) >*/
i__3 = *p - l;
caxpy_(&i__3, &t, &v[lp1 + l * v_dim1], &c__1, &v[lp1 + j *
v_dim1], &c__1);
/*< 310 continue >*/
/* L310: */
}
/*< 320 continue >*/
L320:
/*< do 330 i = 1, p >*/
i__2 = *p;
for (i__ = 1; i__ <= i__2; ++i__) {
/*< v(i,l) = (0.0e0,0.0e0) >*/
i__3 = i__ + l * v_dim1;
v[i__3].r = (float)0., v[i__3].i = (float)0.;
/*< 330 continue >*/
/* L330: */
}
/*< v(l,l) = (1.0e0,0.0e0) >*/
i__2 = l + l * v_dim1;
v[i__2].r = (float)1., v[i__2].i = (float)0.;
/*< 340 continue >*/
/* L340: */
}
/*< 350 continue >*/
L350:
/* transform s and e so that they are real. */
/*< do 380 i = 1, m >*/
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< if (cabs1(s(i)) .eq. 0.0e0) go to 360 >*/
i__2 = i__;
if ((r__1 = s[i__2].r, dabs(r__1)) + (r__2 = r_imag(&s[i__]), dabs(
r__2)) == (float)0.) {
goto L360;
}
/*< t = cmplx(cabs(s(i)),0.0e0) >*/
r__1 = c_abs(&s[i__]);
q__1.r = r__1, q__1.i = (float)0.;
t.r = q__1.r, t.i = q__1.i;
/*< r = s(i)/t >*/
c_div(&q__1, &s[i__], &t);
r__.r = q__1.r, r__.i = q__1.i;
/*< s(i) = t >*/
i__2 = i__;
s[i__2].r = t.r, s[i__2].i = t.i;
/*< if (i .lt. m) e(i) = e(i)/r >*/
if (i__ < m) {
i__2 = i__;
c_div(&q__1, &e[i__], &r__);
e[i__2].r = q__1.r, e[i__2].i = q__1.i;
}
/*< if (wantu) call cscal(n,r,u(1,i),1) >*/
if (wantu) {
cscal_(n, &r__, &u[i__ * u_dim1 + 1], &c__1);
}
/*< 360 continue >*/
L360:
/* ...exit */
/*< if (i .eq. m) go to 390 >*/
if (i__ == m) {
goto L390;
}
/*< if (cabs1(e(i)) .eq. 0.0e0) go to 370 >*/
i__2 = i__;
if ((r__1 = e[i__2].r, dabs(r__1)) + (r__2 = r_imag(&e[i__]), dabs(
r__2)) == (float)0.) {
goto L370;
}
/*< t = cmplx(cabs(e(i)),0.0e0) >*/
r__1 = c_abs(&e[i__]);
q__1.r = r__1, q__1.i = (float)0.;
t.r = q__1.r, t.i = q__1.i;
/*< r = t/e(i) >*/
c_div(&q__1, &t, &e[i__]);
r__.r = q__1.r, r__.i = q__1.i;
/*< e(i) = t >*/
i__2 = i__;
e[i__2].r = t.r, e[i__2].i = t.i;
/*< s(i+1) = s(i+1)*r >*/
i__2 = i__ + 1;
i__3 = i__ + 1;
q__1.r = s[i__3].r * r__.r - s[i__3].i * r__.i, q__1.i = s[i__3].r *
r__.i + s[i__3].i * r__.r;
s[i__2].r = q__1.r, s[i__2].i = q__1.i;
/*< if (wantv) call cscal(p,r,v(1,i+1),1) >*/
if (wantv) {
cscal_(p, &r__, &v[(i__ + 1) * v_dim1 + 1], &c__1);
}
/*< 370 continue >*/
L370:
/*< 380 continue >*/
/* L380: */
;
}
/*< 390 continue >*/
L390:
/* main iteration loop for the singular values. */
/*< mm = m >*/
mm = m;
/*< iter = 0 >*/
iter = 0;
/*< 400 continue >*/
L400:
/* quit if all the singular values have been found. */
/* ...exit */
/*< if (m .eq. 0) go to 660 >*/
if (m == 0) {
goto L660;
}
/* if too many iterations have been performed, set */
/* flag and return. */
/*< if (iter .lt. maxit) go to 410 >*/
if (iter < maxit) {
goto L410;
}
/*< info = m >*/
*info = m;
/* ......exit */
/*< go to 660 >*/
goto L660;
/*< 410 continue >*/
L410:
/* this section of the program inspects for */
/* negligible elements in the s and e arrays. on */
/* completion the variables kase and l are set as follows. */
/* kase = 1 if s(m) and e(l-1) are negligible and l.lt.m */
/* kase = 2 if s(l) is negligible and l.lt.m */
/* kase = 3 if e(l-1) is negligible, l.lt.m, and */
/* s(l), ..., s(m) are not negligible (qr step). */
/* kase = 4 if e(m-1) is negligible (convergence). */
/*< do 430 ll = 1, m >*/
i__1 = m;
for (ll = 1; ll <= i__1; ++ll) {
/*< l = m - ll >*/
l = m - ll;
/* ...exit */
/*< if (l .eq. 0) go to 440 >*/
if (l == 0) {
goto L440;
}
/*< test = cabs(s(l)) + cabs(s(l+1)) >*/
test = c_abs(&s[l]) + c_abs(&s[l + 1]);
/*< ztest = test + cabs(e(l)) >*/
ztest = test + c_abs(&e[l]);
/*< if (ztest .ne. test) go to 420 >*/
if (ztest != test) {
goto L420;
}
/*< e(l) = (0.0e0,0.0e0) >*/
i__2 = l;
e[i__2].r = (float)0., e[i__2].i = (float)0.;
/* ......exit */
/*< go to 440 >*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -