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