?? clals0.c
字號:
Step (1L): apply back the Givens rotations performed. */
latime_1.ops += (real) (*nrhs * 12 * *givptr);
i__1 = *givptr;
for (i__ = 1; i__ <= i__1; ++i__) {
csrot_(nrhs, &b_ref(givcol_ref(i__, 2), 1), ldb, &b_ref(
givcol_ref(i__, 1), 1), ldb, &givnum_ref(i__, 2), &
givnum_ref(i__, 1));
/* L10: */
}
/* Step (2L): permute rows of B. */
ccopy_(nrhs, &b_ref(nlp1, 1), ldb, &bx_ref(1, 1), ldbx);
i__1 = n;
for (i__ = 2; i__ <= i__1; ++i__) {
ccopy_(nrhs, &b_ref(perm[i__], 1), ldb, &bx_ref(i__, 1), ldbx);
/* L20: */
}
/* Step (3L): apply the inverse of the left singular vector
matrix to BX. */
if (*k == 1) {
ccopy_(nrhs, &bx[bx_offset], ldbx, &b[b_offset], ldb);
if (z__[1] < 0.f) {
latime_1.ops += (real) (*nrhs * 6);
csscal_(nrhs, &c_b5, &b[b_offset], ldb);
}
} else {
i__1 = *k;
for (j = 1; j <= i__1; ++j) {
diflj = difl[j];
dj = poles_ref(j, 1);
dsigj = -poles_ref(j, 2);
if (j < *k) {
difrj = -difr_ref(j, 1);
dsigjp = -poles_ref(j + 1, 2);
}
if (z__[j] == 0.f || poles_ref(j, 2) == 0.f) {
rwork[j] = 0.f;
} else {
latime_1.ops += 4.f;
rwork[j] = -poles_ref(j, 2) * z__[j] / diflj / (poles_ref(
j, 2) + dj);
}
i__2 = j - 1;
for (i__ = 1; i__ <= i__2; ++i__) {
if (z__[i__] == 0.f || poles_ref(i__, 2) == 0.f) {
rwork[i__] = 0.f;
} else {
latime_1.ops += 6.f;
rwork[i__] = poles_ref(i__, 2) * z__[i__] / (slamc3_(&
poles_ref(i__, 2), &dsigj) - diflj) / (
poles_ref(i__, 2) + dj);
}
/* L30: */
}
i__2 = *k;
for (i__ = j + 1; i__ <= i__2; ++i__) {
if (z__[i__] == 0.f || poles_ref(i__, 2) == 0.f) {
rwork[i__] = 0.f;
} else {
latime_1.ops += 6.f;
rwork[i__] = poles_ref(i__, 2) * z__[i__] / (slamc3_(&
poles_ref(i__, 2), &dsigjp) + difrj) / (
poles_ref(i__, 2) + dj);
}
/* L40: */
}
latime_1.ops += (real) (*k << 1);
rwork[1] = -1.f;
temp = snrm2_(k, &rwork[1], &c__1);
/* Since B and BX are complex, the following call to SGEMV
is performed in two steps (real and imaginary parts).
CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
$ B( J, 1 ), LDB ) */
i__ = *k + (*nrhs << 1);
i__2 = *nrhs;
for (jcol = 1; jcol <= i__2; ++jcol) {
i__3 = *k;
for (jrow = 1; jrow <= i__3; ++jrow) {
++i__;
i__4 = bx_subscr(jrow, jcol);
rwork[i__] = bx[i__4].r;
/* L50: */
}
/* L60: */
}
latime_1.ops += sopbl2_("SGEMV ", k, nrhs, &c__0, &c__0);
sgemv_("T", k, nrhs, &c_b16, &rwork[*k + 1 + (*nrhs << 1)], k,
&rwork[1], &c__1, &c_b18, &rwork[*k + 1], &c__1);
i__ = *k + (*nrhs << 1);
i__2 = *nrhs;
for (jcol = 1; jcol <= i__2; ++jcol) {
i__3 = *k;
for (jrow = 1; jrow <= i__3; ++jrow) {
++i__;
rwork[i__] = r_imag(&bx_ref(jrow, jcol));
/* L70: */
}
/* L80: */
}
latime_1.ops += sopbl2_("SGEMV ", k, nrhs, &c__0, &c__0);
sgemv_("T", k, nrhs, &c_b16, &rwork[*k + 1 + (*nrhs << 1)], k,
&rwork[1], &c__1, &c_b18, &rwork[*k + 1 + *nrhs], &
c__1);
i__2 = *nrhs;
for (jcol = 1; jcol <= i__2; ++jcol) {
i__3 = b_subscr(j, jcol);
i__4 = jcol + *k;
i__5 = jcol + *k + *nrhs;
q__1.r = rwork[i__4], q__1.i = rwork[i__5];
b[i__3].r = q__1.r, b[i__3].i = q__1.i;
/* L90: */
}
latime_1.ops += (real) (*nrhs * 6);
clascl_("G", &c__0, &c__0, &temp, &c_b16, &c__1, nrhs, &b_ref(
j, 1), ldb, info);
/* L100: */
}
}
/* Move the deflated rows of BX to B also. */
if (*k < max(m,n)) {
i__1 = n - *k;
clacpy_("A", &i__1, nrhs, &bx_ref(*k + 1, 1), ldbx, &b_ref(*k + 1,
1), ldb);
}
} else {
/* Apply back the right orthogonal transformations.
Step (1R): apply back the new right singular vector matrix
to B. */
if (*k == 1) {
ccopy_(nrhs, &b[b_offset], ldb, &bx[bx_offset], ldbx);
} else {
i__1 = *k;
for (j = 1; j <= i__1; ++j) {
dsigj = poles_ref(j, 2);
if (z__[j] == 0.f) {
rwork[j] = 0.f;
} else {
latime_1.ops += 4.f;
rwork[j] = -z__[j] / difl[j] / (dsigj + poles_ref(j, 1)) /
difr_ref(j, 2);
}
i__2 = j - 1;
for (i__ = 1; i__ <= i__2; ++i__) {
if (z__[j] == 0.f) {
rwork[i__] = 0.f;
} else {
latime_1.ops += 6.f;
r__1 = -poles_ref(i__ + 1, 2);
rwork[i__] = z__[j] / (slamc3_(&dsigj, &r__1) -
difr_ref(i__, 1)) / (dsigj + poles_ref(i__, 1)
) / difr_ref(i__, 2);
}
/* L110: */
}
i__2 = *k;
for (i__ = j + 1; i__ <= i__2; ++i__) {
if (z__[j] == 0.f) {
rwork[i__] = 0.f;
} else {
latime_1.ops += 6.f;
r__1 = -poles_ref(i__, 2);
rwork[i__] = z__[j] / (slamc3_(&dsigj, &r__1) - difl[
i__]) / (dsigj + poles_ref(i__, 1)) /
difr_ref(i__, 2);
}
/* L120: */
}
/* Since B and BX are complex, the following call to SGEMV
is performed in two steps (real and imaginary parts).
CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
$ BX( J, 1 ), LDBX ) */
i__ = *k + (*nrhs << 1);
i__2 = *nrhs;
for (jcol = 1; jcol <= i__2; ++jcol) {
i__3 = *k;
for (jrow = 1; jrow <= i__3; ++jrow) {
++i__;
i__4 = b_subscr(jrow, jcol);
rwork[i__] = b[i__4].r;
/* L130: */
}
/* L140: */
}
latime_1.ops += sopbl2_("SGEMV ", k, nrhs, &c__0, &c__0);
sgemv_("T", k, nrhs, &c_b16, &rwork[*k + 1 + (*nrhs << 1)], k,
&rwork[1], &c__1, &c_b18, &rwork[*k + 1], &c__1);
i__ = *k + (*nrhs << 1);
i__2 = *nrhs;
for (jcol = 1; jcol <= i__2; ++jcol) {
i__3 = *k;
for (jrow = 1; jrow <= i__3; ++jrow) {
++i__;
rwork[i__] = r_imag(&b_ref(jrow, jcol));
/* L150: */
}
/* L160: */
}
latime_1.ops += sopbl2_("SGEMV ", k, nrhs, &c__0, &c__0);
sgemv_("T", k, nrhs, &c_b16, &rwork[*k + 1 + (*nrhs << 1)], k,
&rwork[1], &c__1, &c_b18, &rwork[*k + 1 + *nrhs], &
c__1);
i__2 = *nrhs;
for (jcol = 1; jcol <= i__2; ++jcol) {
i__3 = bx_subscr(j, jcol);
i__4 = jcol + *k;
i__5 = jcol + *k + *nrhs;
q__1.r = rwork[i__4], q__1.i = rwork[i__5];
bx[i__3].r = q__1.r, bx[i__3].i = q__1.i;
/* L170: */
}
/* L180: */
}
}
/* Step (2R): if SQRE = 1, apply back the rotation that is
related to the right null space of the subproblem. */
if (*sqre == 1) {
latime_1.ops += (real) (*nrhs * 12);
ccopy_(nrhs, &b_ref(m, 1), ldb, &bx_ref(m, 1), ldbx);
csrot_(nrhs, &bx_ref(1, 1), ldbx, &bx_ref(m, 1), ldbx, c__, s);
}
if (*k < max(m,n)) {
i__1 = n - *k;
clacpy_("A", &i__1, nrhs, &b_ref(*k + 1, 1), ldb, &bx_ref(*k + 1,
1), ldbx);
}
/* Step (3R): permute rows of B. */
ccopy_(nrhs, &bx_ref(1, 1), ldbx, &b_ref(nlp1, 1), ldb);
if (*sqre == 1) {
ccopy_(nrhs, &bx_ref(m, 1), ldbx, &b_ref(m, 1), ldb);
}
i__1 = n;
for (i__ = 2; i__ <= i__1; ++i__) {
ccopy_(nrhs, &bx_ref(i__, 1), ldbx, &b_ref(perm[i__], 1), ldb);
/* L190: */
}
/* Step (4R): apply back the Givens rotations performed. */
latime_1.ops += (real) (*nrhs * 12 * *givptr);
for (i__ = *givptr; i__ >= 1; --i__) {
r__1 = -givnum_ref(i__, 1);
csrot_(nrhs, &b_ref(givcol_ref(i__, 2), 1), ldb, &b_ref(
givcol_ref(i__, 1), 1), ldb, &givnum_ref(i__, 2), &r__1);
/* L200: */
}
}
return 0;
/* End of CLALS0 */
} /* clals0_ */
#undef givnum_ref
#undef givcol_ref
#undef bx_ref
#undef bx_subscr
#undef poles_ref
#undef b_ref
#undef b_subscr
#undef difr_ref
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -