📄 stgsja.c
字号:
} else if (! (initq || wantq || lsame_(jobq, "N", (ftnlen)1, (ftnlen)1)))
{
/*< INFO = -3 >*/
*info = -3;
/*< ELSE IF( M.LT.0 ) THEN >*/
} else if (*m < 0) {
/*< INFO = -4 >*/
*info = -4;
/*< ELSE IF( P.LT.0 ) THEN >*/
} else if (*p < 0) {
/*< INFO = -5 >*/
*info = -5;
/*< ELSE IF( N.LT.0 ) THEN >*/
} else if (*n < 0) {
/*< INFO = -6 >*/
*info = -6;
/*< ELSE IF( LDA.LT.MAX( 1, M ) ) THEN >*/
} else if (*lda < max(1,*m)) {
/*< INFO = -10 >*/
*info = -10;
/*< ELSE IF( LDB.LT.MAX( 1, P ) ) THEN >*/
} else if (*ldb < max(1,*p)) {
/*< INFO = -12 >*/
*info = -12;
/*< ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN >*/
} else if (*ldu < 1 || (wantu && *ldu < *m)) {
/*< INFO = -18 >*/
*info = -18;
/*< ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN >*/
} else if (*ldv < 1 || (wantv && *ldv < *p)) {
/*< INFO = -20 >*/
*info = -20;
/*< ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN >*/
} else if (*ldq < 1 || (wantq && *ldq < *n)) {
/*< INFO = -22 >*/
*info = -22;
/*< END IF >*/
}
/*< IF( INFO.NE.0 ) THEN >*/
if (*info != 0) {
/*< CALL XERBLA( 'STGSJA', -INFO ) >*/
i__1 = -(*info);
xerbla_("STGSJA", &i__1, (ftnlen)6);
/*< RETURN >*/
return 0;
/*< END IF >*/
}
/* Initialize U, V and Q, if necessary */
/*< >*/
if (initu) {
slaset_("Full", m, m, &c_b13, &c_b14, &u[u_offset], ldu, (ftnlen)4);
}
/*< >*/
if (initv) {
slaset_("Full", p, p, &c_b13, &c_b14, &v[v_offset], ldv, (ftnlen)4);
}
/*< >*/
if (initq) {
slaset_("Full", n, n, &c_b13, &c_b14, &q[q_offset], ldq, (ftnlen)4);
}
/* Loop until convergence */
/*< UPPER = .FALSE. >*/
upper = FALSE_;
/*< DO 40 KCYCLE = 1, MAXIT >*/
for (kcycle = 1; kcycle <= 40; ++kcycle) {
/*< UPPER = .NOT.UPPER >*/
upper = ! upper;
/*< DO 20 I = 1, L - 1 >*/
i__1 = *l - 1;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< DO 10 J = I + 1, L >*/
i__2 = *l;
for (j = i__ + 1; j <= i__2; ++j) {
/*< A1 = ZERO >*/
a1 = (float)0.;
/*< A2 = ZERO >*/
a2 = (float)0.;
/*< A3 = ZERO >*/
a3 = (float)0.;
/*< >*/
if (*k + i__ <= *m) {
a1 = a[*k + i__ + (*n - *l + i__) * a_dim1];
}
/*< >*/
if (*k + j <= *m) {
a3 = a[*k + j + (*n - *l + j) * a_dim1];
}
/*< B1 = B( I, N-L+I ) >*/
b1 = b[i__ + (*n - *l + i__) * b_dim1];
/*< B3 = B( J, N-L+J ) >*/
b3 = b[j + (*n - *l + j) * b_dim1];
/*< IF( UPPER ) THEN >*/
if (upper) {
/*< >*/
if (*k + i__ <= *m) {
a2 = a[*k + i__ + (*n - *l + j) * a_dim1];
}
/*< B2 = B( I, N-L+J ) >*/
b2 = b[i__ + (*n - *l + j) * b_dim1];
/*< ELSE >*/
} else {
/*< >*/
if (*k + j <= *m) {
a2 = a[*k + j + (*n - *l + i__) * a_dim1];
}
/*< B2 = B( J, N-L+I ) >*/
b2 = b[j + (*n - *l + i__) * b_dim1];
/*< END IF >*/
}
/*< >*/
slags2_(&upper, &a1, &a2, &a3, &b1, &b2, &b3, &csu, &snu, &
csv, &snv, &csq, &snq);
/* Update (K+I)-th and (K+J)-th rows of matrix A: U'*A */
/*< >*/
if (*k + j <= *m) {
srot_(l, &a[*k + j + (*n - *l + 1) * a_dim1], lda, &a[*k
+ i__ + (*n - *l + 1) * a_dim1], lda, &csu, &snu);
}
/* Update I-th and J-th rows of matrix B: V'*B */
/*< >*/
srot_(l, &b[j + (*n - *l + 1) * b_dim1], ldb, &b[i__ + (*n - *
l + 1) * b_dim1], ldb, &csv, &snv);
/* Update (N-L+I)-th and (N-L+J)-th columns of matrices */
/* A and B: A*Q and B*Q */
/*< >*/
/* Computing MIN */
i__4 = *k + *l;
i__3 = min(i__4,*m);
srot_(&i__3, &a[(*n - *l + j) * a_dim1 + 1], &c__1, &a[(*n - *
l + i__) * a_dim1 + 1], &c__1, &csq, &snq);
/*< >*/
srot_(l, &b[(*n - *l + j) * b_dim1 + 1], &c__1, &b[(*n - *l +
i__) * b_dim1 + 1], &c__1, &csq, &snq);
/*< IF( UPPER ) THEN >*/
if (upper) {
/*< >*/
if (*k + i__ <= *m) {
a[*k + i__ + (*n - *l + j) * a_dim1] = (float)0.;
}
/*< B( I, N-L+J ) = ZERO >*/
b[i__ + (*n - *l + j) * b_dim1] = (float)0.;
/*< ELSE >*/
} else {
/*< >*/
if (*k + j <= *m) {
a[*k + j + (*n - *l + i__) * a_dim1] = (float)0.;
}
/*< B( J, N-L+I ) = ZERO >*/
b[j + (*n - *l + i__) * b_dim1] = (float)0.;
/*< END IF >*/
}
/* Update orthogonal matrices U, V, Q, if desired. */
/*< >*/
if (wantu && *k + j <= *m) {
srot_(m, &u[(*k + j) * u_dim1 + 1], &c__1, &u[(*k + i__) *
u_dim1 + 1], &c__1, &csu, &snu);
}
/*< >*/
if (wantv) {
srot_(p, &v[j * v_dim1 + 1], &c__1, &v[i__ * v_dim1 + 1],
&c__1, &csv, &snv);
}
/*< >*/
if (wantq) {
srot_(n, &q[(*n - *l + j) * q_dim1 + 1], &c__1, &q[(*n - *
l + i__) * q_dim1 + 1], &c__1, &csq, &snq);
}
/*< 10 CONTINUE >*/
/* L10: */
}
/*< 20 CONTINUE >*/
/* L20: */
}
/*< IF( .NOT.UPPER ) THEN >*/
if (! upper) {
/* The matrices A13 and B13 were lower triangular at the start */
/* of the cycle, and are now upper triangular. */
/* Convergence test: test the parallelism of the corresponding */
/* rows of A and B. */
/*< ERROR = ZERO >*/
error = (float)0.;
/*< DO 30 I = 1, MIN( L, M-K ) >*/
/* Computing MIN */
i__2 = *l, i__3 = *m - *k;
i__1 = min(i__2,i__3);
for (i__ = 1; i__ <= i__1; ++i__) {
/*< CALL SCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 ) >*/
i__2 = *l - i__ + 1;
scopy_(&i__2, &a[*k + i__ + (*n - *l + i__) * a_dim1], lda, &
work[1], &c__1);
/*< CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 ) >*/
i__2 = *l - i__ + 1;
scopy_(&i__2, &b[i__ + (*n - *l + i__) * b_dim1], ldb, &work[*
l + 1], &c__1);
/*< CALL SLAPLL( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN ) >*/
i__2 = *l - i__ + 1;
slapll_(&i__2, &work[1], &c__1, &work[*l + 1], &c__1, &ssmin);
/*< ERROR = MAX( ERROR, SSMIN ) >*/
error = dmax(error,ssmin);
/*< 30 CONTINUE >*/
/* L30: */
}
/*< >*/
if (dabs(error) <= dmin(*tola,*tolb)) {
goto L50;
}
/*< END IF >*/
}
/* End of cycle loop */
/*< 40 CONTINUE >*/
/* L40: */
}
/* The algorithm has not converged after MAXIT cycles. */
/*< INFO = 1 >*/
*info = 1;
/*< GO TO 100 >*/
goto L100;
/*< 50 CONTINUE >*/
L50:
/* If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged. */
/* Compute the generalized singular value pairs (ALPHA, BETA), and */
/* set the triangular matrix R to array A. */
/*< DO 60 I = 1, K >*/
i__1 = *k;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< ALPHA( I ) = ONE >*/
alpha[i__] = (float)1.;
/*< BETA( I ) = ZERO >*/
beta[i__] = (float)0.;
/*< 60 CONTINUE >*/
/* L60: */
}
/*< DO 70 I = 1, MIN( L, M-K ) >*/
/* Computing MIN */
i__2 = *l, i__3 = *m - *k;
i__1 = min(i__2,i__3);
for (i__ = 1; i__ <= i__1; ++i__) {
/*< A1 = A( K+I, N-L+I ) >*/
a1 = a[*k + i__ + (*n - *l + i__) * a_dim1];
/*< B1 = B( I, N-L+I ) >*/
b1 = b[i__ + (*n - *l + i__) * b_dim1];
/*< IF( A1.NE.ZERO ) THEN >*/
if (a1 != (float)0.) {
/*< GAMMA = B1 / A1 >*/
gamma = b1 / a1;
/* change sign if necessary */
/*< IF( GAMMA.LT.ZERO ) THEN >*/
if (gamma < (float)0.) {
/*< CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) >*/
i__2 = *l - i__ + 1;
sscal_(&i__2, &c_b43, &b[i__ + (*n - *l + i__) * b_dim1], ldb)
;
/*< >*/
if (wantv) {
sscal_(p, &c_b43, &v[i__ * v_dim1 + 1], &c__1);
}
/*< END IF >*/
}
/*< >*/
r__1 = dabs(gamma);
slartg_(&r__1, &c_b14, &beta[*k + i__], &alpha[*k + i__], &rwk);
/*< IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN >*/
if (alpha[*k + i__] >= beta[*k + i__]) {
/*< >*/
i__2 = *l - i__ + 1;
r__1 = (float)1. / alpha[*k + i__];
sscal_(&i__2, &r__1, &a[*k + i__ + (*n - *l + i__) * a_dim1],
lda);
/*< ELSE >*/
} else {
/*< >*/
i__2 = *l - i__ + 1;
r__1 = (float)1. / beta[*k + i__];
sscal_(&i__2, &r__1, &b[i__ + (*n - *l + i__) * b_dim1], ldb);
/*< >*/
i__2 = *l - i__ + 1;
scopy_(&i__2, &b[i__ + (*n - *l + i__) * b_dim1], ldb, &a[*k
+ i__ + (*n - *l + i__) * a_dim1], lda);
/*< END IF >*/
}
/*< ELSE >*/
} else {
/*< ALPHA( K+I ) = ZERO >*/
alpha[*k + i__] = (float)0.;
/*< BETA( K+I ) = ONE >*/
beta[*k + i__] = (float)1.;
/*< >*/
i__2 = *l - i__ + 1;
scopy_(&i__2, &b[i__ + (*n - *l + i__) * b_dim1], ldb, &a[*k +
i__ + (*n - *l + i__) * a_dim1], lda);
/*< END IF >*/
}
/*< 70 CONTINUE >*/
/* L70: */
}
/* Post-assignment */
/*< DO 80 I = M + 1, K + L >*/
i__1 = *k + *l;
for (i__ = *m + 1; i__ <= i__1; ++i__) {
/*< ALPHA( I ) = ZERO >*/
alpha[i__] = (float)0.;
/*< BETA( I ) = ONE >*/
beta[i__] = (float)1.;
/*< 80 CONTINUE >*/
/* L80: */
}
/*< IF( K+L.LT.N ) THEN >*/
if (*k + *l < *n) {
/*< DO 90 I = K + L + 1, N >*/
i__1 = *n;
for (i__ = *k + *l + 1; i__ <= i__1; ++i__) {
/*< ALPHA( I ) = ZERO >*/
alpha[i__] = (float)0.;
/*< BETA( I ) = ZERO >*/
beta[i__] = (float)0.;
/*< 90 CONTINUE >*/
/* L90: */
}
/*< END IF >*/
}
/*< 100 CONTINUE >*/
L100:
/*< NCYCLE = KCYCLE >*/
*ncycle = kcycle;
/*< RETURN >*/
return 0;
/* End of STGSJA */
/*< END >*/
} /* stgsja_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -