dtgsja.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 505 行 · 第 1/2 页
C
505 行
/* */
/* MAXIT INTEGER */
/* MAXIT specifies the total loops that the iterative procedure */
/* may take. If after MAXIT cycles, the routine fails to */
/* converge, we return INFO = 1. */
/* */
/* Further Details */
/* =============== */
/* */
/* DTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce */
/* min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L */
/* matrix B13 to the form: */
/* */
/* U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1, */
/* */
/* where U1, V1 and Q1 are orthogonal matrix, and Z' is the transpose */
/* of Z. C1 and S1 are diagonal matrices satisfying */
/* */
/* C1**2 + S1**2 = I, */
/* */
/* and R1 is an L-by-L nonsingular upper triangular matrix. */
/* */
/* ===================================================================== */
/* Decode and test the input parameters */
initu = lsame_(jobu, "I");
wantu = initu || lsame_(jobu, "U");
initv = lsame_(jobv, "I");
wantv = initv || lsame_(jobv, "V");
initq = lsame_(jobq, "I");
wantq = initq || lsame_(jobq, "Q");
*info = 0;
if (! (initu || wantu || lsame_(jobu, "N"))) {
*info = -1;
} else if (! (initv || wantv || lsame_(jobv, "N"))) {
*info = -2;
} else if (! (initq || wantq || lsame_(jobq, "N"))) {
*info = -3;
} else if (*m < 0) {
*info = -4;
} else if (*p < 0) {
*info = -5;
} else if (*n < 0) {
*info = -6;
} else if (*lda < max(1,*m)) {
*info = -10;
} else if (*ldb < max(1,*p)) {
*info = -12;
} else if (*ldu < 1 || (wantu && *ldu < *m)) {
*info = -18;
} else if (*ldv < 1 || (wantv && *ldv < *p)) {
*info = -20;
} else if (*ldq < 1 || (wantq && *ldq < *n)) {
*info = -22;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DTGSJA", &i__1);
return;
}
/* Initialize U, V and Q, if necessary */
if (initu) {
dlaset_("Full", m, m, &c_b13, &c_b14, u, ldu);
}
if (initv) {
dlaset_("Full", p, p, &c_b13, &c_b14, v, ldv);
}
if (initq) {
dlaset_("Full", n, n, &c_b13, &c_b14, q, ldq);
}
/* Loop until convergence */
upper = FALSE_;
for (kcycle = 1; kcycle <= 40; ++kcycle) {
upper = ! upper;
for (i = 0; i < *l; ++i) {
for (j = i + 1; j < *l; ++j) {
a1 = a2 = a3 = 0.;
if (*k + i < *m) {
a1 = a[*k + i + (*n - *l + i) * *lda];
}
if (*k + j < *m) {
a3 = a[*k + j + (*n - *l + j) * *lda];
}
b1 = b[i + (*n - *l + i) * *ldb];
b3 = b[j + (*n - *l + j) * *ldb];
if (upper) {
if (*k + i < *m) {
a2 = a[*k + i + (*n - *l + j) * *lda];
}
b2 = b[i + (*n - *l + j) * *ldb];
} else {
if (*k + j < *m) {
a2 = a[*k + j + (*n - *l + i) * *lda];
}
b2 = b[j + (*n - *l + i) * *ldb];
}
dlags2_(&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) {
drot_(l, &a[*k + j + (*n - *l) * *lda], lda, &a[*k + i + (*n - *l) * *lda], lda, &csu, &snu);
}
/* Update I-th and J-th rows of matrix B: V'*B */
drot_(l, &b[j + (*n - *l) * *ldb], ldb, &b[i + (*n - *l) * *ldb], ldb, &csv, &snv);
/* Update (N-L+I)-th and (N-L+J)-th columns of matrices */
/* A and B: A*Q and B*Q */
i__1 = min(*k + *l,*m);
drot_(&i__1, &a[(*n - *l + j) * *lda], &c__1, &a[(*n - *l + i) * *lda], &c__1, &csq, &snq);
drot_(l, &b[(*n - *l + j) * *ldb], &c__1, &b[(*n - *l + i) * *ldb], &c__1, &csq, &snq);
if (upper) {
if (*k + i < *m) {
a[*k + i + (*n - *l + j) * *lda] = 0.;
}
b[i + (*n - *l + j) * *ldb] = 0.;
} else {
if (*k + j < *m) {
a[*k + j + (*n - *l + i) * *lda] = 0.;
}
b[j + (*n - *l + i) * *ldb] = 0.;
}
/* Update orthogonal matrices U, V, Q, if desired. */
if (wantu && *k + j < *m) {
drot_(m, &u[(*k + j) * *ldu], &c__1, &u[(*k + i) * *ldu], &c__1, &csu, &snu);
}
if (wantv) {
drot_(p, &v[j * *ldv], &c__1, &v[i * *ldv], &c__1, &csv, &snv);
}
if (wantq) {
drot_(n, &q[(*n - *l + j) * *ldq], &c__1, &q[(*n - *l + i) * *ldq], &c__1, &csq, &snq);
}
}
}
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 = 0.;
for (i = 0; i < min(*l,*m - *k); ++i) {
i__1 = *l - i;
dcopy_(&i__1, &a[*k + i + (*n - *l + i) * *lda], lda, work, &c__1);
dcopy_(&i__1, &b[i + (*n - *l + i) * *ldb], ldb, &work[*l], &c__1);
dlapll_(&i__1, work, &c__1, &work[*l], &c__1, &ssmin);
error = max(error,ssmin);
}
if (abs(error) <= (doublereal) (*n) * min(*tola,*tolb)) {
goto L50;
}
}
/* End of cycle loop */
}
/* The algorithm has not converged after MAXIT cycles. */
*info = 1;
*ncycle = kcycle;
return;
L50:
/* If ERROR <= N*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. */
for (i = 0; i < *k; ++i) {
alpha[i] = 1.;
beta[i] = 0.;
}
for (i = 0; i < min(*l,*m - *k); ++i) {
a1 = a[*k + i + (*n - *l + i) * *lda];
b1 = b[i + (*n - *l + i) * *ldb];
if (a1 != 0.) {
gamma = b1 / a1;
/* change sign if necessary */
if (gamma < 0.) {
i__1 = *l - i;
dscal_(&i__1, &c_b43, &b[i + (*n - *l + i) * *ldb], ldb);
if (wantv) {
dscal_(p, &c_b43, &v[i * *ldv], &c__1);
}
}
d__1 = abs(gamma);
dlartg_(&d__1, &c_b14, &beta[*k + i], &alpha[*k + i], &rwk);
if (alpha[*k + i] >= beta[*k + i]) {
i__1 = *l - i;
d__1 = 1. / alpha[*k + i];
dscal_(&i__1, &d__1, &a[*k + i + (*n - *l + i) * *lda], lda);
} else {
i__1 = *l - i;
d__1 = 1. / beta[*k + i];
dscal_(&i__1, &d__1, &b[i + (*n - *l + i) * *ldb], ldb);
dcopy_(&i__1, &b[i + (*n - *l + i) * *ldb], ldb, &a[*k + i + (*n - *l + i) * *lda], lda);
}
} else {
alpha[*k + i] = 0.;
beta[*k + i] = 1.;
i__1 = *l - i;
dcopy_(&i__1, &b[i + (*n - *l + i) * *ldb], ldb, &a[*k + i + (*n - *l + i) * *lda], lda);
}
}
/* Post-assignment */
for (i = *m; i < *k + *l; ++i) {
alpha[i] = 0.;
beta[i] = 1.;
}
if (*k + *l < *n) {
for (i = *k + *l; i < *n; ++i) {
alpha[i] = 0.;
beta[i] = 0.;
}
}
*ncycle = kcycle;
return;
} /* dtgsja_ */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?