dtgsen.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 752 行 · 第 1/3 页
C
752 行
/* decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the */
/* reordered generalized real Schur form of (C, D) is given by */
/* */
/* (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)', */
/* */
/* and the first n1 columns of Q*U and Z*W span the corresponding */
/* deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.). */
/* */
/* Note that if the selected eigenvalue is sufficiently ill-conditioned, */
/* then its value may differ significantly from its value before */
/* reordering. */
/* */
/* The reciprocal condition numbers of the left and right eigenspaces */
/* spanned by the first n1 columns of U and W (or Q*U and Z*W) may */
/* be returned in DIF(1:2), corresponding to Difu and Difl, resp. */
/* */
/* The Difu and Difl are defined as: */
/* */
/* Difu[(A11, B11), (A22, B22)] = sigma-min( Zu ) */
/* and */
/* Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)], */
/* */
/* where sigma-min(Zu) is the smallest singular value of the */
/* (2*n1*n2)-by-(2*n1*n2) matrix */
/* */
/* Zu = [ kron(In2, A11) -kron(A22', In1) ] */
/* [ kron(In2, B11) -kron(B22', In1) ]. */
/* */
/* Here, Inx is the identity matrix of size nx and A22' is the */
/* transpose of A22. kron(X, Y) is the Kronecker product between */
/* the matrices X and Y. */
/* */
/* When DIF(2) is small, small changes in (A, B) can cause large changes */
/* in the deflating subspace. An approximate (asymptotic) bound on the */
/* maximum angular error in the computed deflating subspaces is */
/* */
/* EPS * norm((A, B)) / DIF(2), */
/* */
/* where EPS is the machine precision. */
/* */
/* The reciprocal norm of the projectors on the left and right */
/* eigenspaces associated with (A11, B11) may be returned in PL and PR. */
/* They are computed as follows. First we compute L and R so that */
/* P*(A, B)*Q is block diagonal, where */
/* */
/* P = ( I -L ) n1 Q = ( I R ) n1 */
/* ( 0 I ) n2 and ( 0 I ) n2 */
/* n1 n2 n1 n2 */
/* */
/* and (L, R) is the solution to the generalized Sylvester equation */
/* */
/* A11*R - L*A22 = -A12 */
/* B11*R - L*B22 = -B12 */
/* */
/* Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2). */
/* An approximate (asymptotic) bound on the average absolute error of */
/* the selected eigenvalues is */
/* */
/* EPS * norm((A, B)) / PL. */
/* */
/* There are also global error bounds which valid for perturbations up */
/* to a certain restriction: A lower bound (x) on the smallest */
/* F-norm(E,F) for which an eigenvalue of (A11, B11) may move and */
/* coalesce with an eigenvalue of (A22, B22) under perturbation (E,F), */
/* (i.e. (A + E, B + F), is */
/* */
/* x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)). */
/* */
/* An approximate bound on x can be computed from DIF(1:2), PL and PR. */
/* */
/* If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed */
/* (L', R') and unperturbed (L, R) left and right deflating subspaces */
/* associated with the selected cluster in the (1,1)-blocks can be */
/* bounded as */
/* */
/* max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2)) */
/* max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2)) */
/* */
/* See LAPACK User's Guide section 4.11 or the following references */
/* for more information. */
/* */
/* Note that if the default method for computing the Frobenius-norm- */
/* based estimate DIF is not wanted (see DLATDF), then the parameter */
/* IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF */
/* (IJOB = 2 will be used)). See DTGSYL for more details. */
/* */
/* Based on contributions by */
/* Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
/* Umea University, S-901 87 Umea, Sweden. */
/* */
/* References */
/* ========== */
/* */
/* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the */
/* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in */
/* M.S. Moonen et al (eds), Linear Algebra for Large Scale and */
/* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. */
/* */
/* [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified */
/* Eigenvalues of a Regular Matrix Pair (A, B) and Condition */
/* Estimation: Theory, Algorithms and Software, */
/* Report UMINF - 94.04, Department of Computing Science, Umea */
/* University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working */
/* Note 87. To appear in Numerical Algorithms, 1996. */
/* */
/* [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software */
/* for Solving the Generalized Sylvester Equation and Estimating the */
/* Separation between Regular Matrix Pairs, Report UMINF - 93.23, */
/* Department of Computing Science, Umea University, S-901 87 Umea, */
/* Sweden, December 1993, Revised April 1994, Also as LAPACK Working */
/* Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, */
/* 1996. */
/* */
/* ===================================================================== */
/* Parameter adjustments */
--select;
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
b_dim1 = *ldb;
b_offset = 1 + b_dim1 * 1;
b -= b_offset;
--alphar;
--alphai;
--beta;
q_dim1 = *ldq;
q_offset = 1 + q_dim1 * 1;
q -= q_offset;
z_dim1 = *ldz;
z_offset = 1 + z_dim1 * 1;
z -= z_offset;
--dif;
--work;
--iwork;
/* Decode and test the input parameters */
*info = 0;
lquery = *lwork == -1 || *liwork == -1;
if (*ijob < 0 || *ijob > 5) {
*info = -1;
} else if (*n < 0) {
*info = -5;
} else if (*lda < max(1,*n)) {
*info = -7;
} else if (*ldb < max(1,*n)) {
*info = -9;
} else if (*ldq < 1 || (*wantq && *ldq < *n)) {
*info = -14;
} else if (*ldz < 1 || (*wantz && *ldz < *n)) {
*info = -16;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DTGSEN", &i__1);
return;
}
/* Get machine constants */
eps = dlamch_("P");
smlnum = dlamch_("S") / eps;
ierr = 0;
wantp = *ijob == 1 || *ijob >= 4;
wantd1 = *ijob == 2 || *ijob == 4;
wantd2 = *ijob == 3 || *ijob == 5;
wantd = wantd1 || wantd2;
/* Set M to the dimension of the specified pair of deflating */
/* subspaces. */
*m = 0;
pair = FALSE_;
for (k = 1; k <= *n; ++k) {
if (pair) {
pair = FALSE_;
} else {
if (k < *n) {
if (a[k + 1 + k * a_dim1] == 0.) {
if (select[k]) {
++(*m);
}
} else {
pair = TRUE_;
if (select[k] || select[k + 1]) {
*m += 2;
}
}
} else {
if (select[*n]) {
++(*m);
}
}
}
}
if (*ijob == 1 || *ijob == 2 || *ijob == 4) {
lwmin = max(max(1, (*n << 2) + 16), (*m << 1) * (*n - *m));
liwmin = max(1, *n + 6);
} else if (*ijob == 3 || *ijob == 5) {
lwmin = max(max(1, (*n << 2) + 16), (*m << 2) * (*n - *m));
liwmin = max(max(1, (*m << 1) * (*n - *m)), *n + 6);
} else {
lwmin = max(1,(*n << 2) + 16);
liwmin = 1;
}
work[1] = (doublereal) lwmin;
iwork[1] = liwmin;
if (*lwork < lwmin && ! lquery) {
*info = -22;
} else if (*liwork < liwmin && ! lquery) {
*info = -24;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DTGSEN", &i__1);
return;
} else if (lquery) {
return;
}
/* Quick return if possible. */
if (*m == *n || *m == 0) {
if (wantp) {
*pl = 1.;
*pr = 1.;
}
if (wantd) {
dscale = 0.;
dsum = 1.;
for (i = 1; i <= *n; ++i) {
dlassq_(n, &a[i * a_dim1 + 1], &c__1, &dscale, &dsum);
dlassq_(n, &b[i * b_dim1 + 1], &c__1, &dscale, &dsum);
}
dif[1] = dscale * sqrt(dsum);
dif[2] = dif[1];
}
goto L60;
}
/* Collect the selected blocks at the top-left corner of (A, B). */
ks = 0;
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?