📄 dtgsy2.c
字号:
/*< INFO = -10 >*/
*info = -10;
/*< ELSE IF( LDD.LT.MAX( 1, M ) ) THEN >*/
} else if (*ldd < max(1,*m)) {
/*< INFO = -12 >*/
*info = -12;
/*< ELSE IF( LDE.LT.MAX( 1, N ) ) THEN >*/
} else if (*lde < max(1,*n)) {
/*< INFO = -14 >*/
*info = -14;
/*< ELSE IF( LDF.LT.MAX( 1, M ) ) THEN >*/
} else if (*ldf < max(1,*m)) {
/*< INFO = -16 >*/
*info = -16;
/*< END IF >*/
}
/*< IF( INFO.NE.0 ) THEN >*/
if (*info != 0) {
/*< CALL XERBLA( 'DTGSY2', -INFO ) >*/
i__1 = -(*info);
xerbla_("DTGSY2", &i__1, (ftnlen)6);
/*< RETURN >*/
return 0;
/*< END IF >*/
}
/* Determine block structure of A */
/*< PQ = 0 >*/
*pq = 0;
/*< P = 0 >*/
p = 0;
/*< I = 1 >*/
i__ = 1;
/*< 10 CONTINUE >*/
L10:
/*< >*/
if (i__ > *m) {
goto L20;
}
/*< P = P + 1 >*/
++p;
/*< IWORK( P ) = I >*/
iwork[p] = i__;
/*< >*/
if (i__ == *m) {
goto L20;
}
/*< IF( A( I+1, I ).NE.ZERO ) THEN >*/
if (a[i__ + 1 + i__ * a_dim1] != 0.) {
/*< I = I + 2 >*/
i__ += 2;
/*< ELSE >*/
} else {
/*< I = I + 1 >*/
++i__;
/*< END IF >*/
}
/*< GO TO 10 >*/
goto L10;
/*< 20 CONTINUE >*/
L20:
/*< IWORK( P+1 ) = M + 1 >*/
iwork[p + 1] = *m + 1;
/* Determine block structure of B */
/*< Q = P + 1 >*/
q = p + 1;
/*< J = 1 >*/
j = 1;
/*< 30 CONTINUE >*/
L30:
/*< >*/
if (j > *n) {
goto L40;
}
/*< Q = Q + 1 >*/
++q;
/*< IWORK( Q ) = J >*/
iwork[q] = j;
/*< >*/
if (j == *n) {
goto L40;
}
/*< IF( B( J+1, J ).NE.ZERO ) THEN >*/
if (b[j + 1 + j * b_dim1] != 0.) {
/*< J = J + 2 >*/
j += 2;
/*< ELSE >*/
} else {
/*< J = J + 1 >*/
++j;
/*< END IF >*/
}
/*< GO TO 30 >*/
goto L30;
/*< 40 CONTINUE >*/
L40:
/*< IWORK( Q+1 ) = N + 1 >*/
iwork[q + 1] = *n + 1;
/*< PQ = P*( Q-P-1 ) >*/
*pq = p * (q - p - 1);
/*< IF( NOTRAN ) THEN >*/
if (notran) {
/* Solve (I, J) - subsystem */
/* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) */
/* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) */
/* for I = P, P - 1, ..., 1; J = 1, 2, ..., Q */
/*< SCALE = ONE >*/
*scale = 1.;
/*< SCALOC = ONE >*/
scaloc = 1.;
/*< DO 120 J = P + 2, Q >*/
i__1 = q;
for (j = p + 2; j <= i__1; ++j) {
/*< JS = IWORK( J ) >*/
js = iwork[j];
/*< JSP1 = JS + 1 >*/
jsp1 = js + 1;
/*< JE = IWORK( J+1 ) - 1 >*/
je = iwork[j + 1] - 1;
/*< NB = JE - JS + 1 >*/
nb = je - js + 1;
/*< DO 110 I = P, 1, -1 >*/
for (i__ = p; i__ >= 1; --i__) {
/*< IS = IWORK( I ) >*/
is = iwork[i__];
/*< ISP1 = IS + 1 >*/
isp1 = is + 1;
/*< IE = IWORK( I+1 ) - 1 >*/
ie = iwork[i__ + 1] - 1;
/*< MB = IE - IS + 1 >*/
mb = ie - is + 1;
/*< ZDIM = MB*NB*2 >*/
zdim = mb * nb << 1;
/*< IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN >*/
if (mb == 1 && nb == 1) {
/* Build a 2-by-2 system Z * x = RHS */
/*< Z( 1, 1 ) = A( IS, IS ) >*/
z__[0] = a[is + is * a_dim1];
/*< Z( 2, 1 ) = D( IS, IS ) >*/
z__[1] = d__[is + is * d_dim1];
/*< Z( 1, 2 ) = -B( JS, JS ) >*/
z__[8] = -b[js + js * b_dim1];
/*< Z( 2, 2 ) = -E( JS, JS ) >*/
z__[9] = -e[js + js * e_dim1];
/* Set up right hand side(s) */
/*< RHS( 1 ) = C( IS, JS ) >*/
rhs[0] = c__[is + js * c_dim1];
/*< RHS( 2 ) = F( IS, JS ) >*/
rhs[1] = f[is + js * f_dim1];
/* Solve Z * x = RHS */
/*< CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) >*/
dgetc2_(&zdim, z__, &c__8, ipiv, jpiv, &ierr);
/*< >*/
if (ierr > 0) {
*info = ierr;
}
/*< IF( IJOB.EQ.0 ) THEN >*/
if (*ijob == 0) {
/*< >*/
dgesc2_(&zdim, z__, &c__8, rhs, ipiv, jpiv, &scaloc);
/*< IF( SCALOC.NE.ONE ) THEN >*/
if (scaloc != 1.) {
/*< DO 50 K = 1, N >*/
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
/*< CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &c__[k * c_dim1 + 1], &
c__1);
/*< CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
/*< 50 CONTINUE >*/
/* L50: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/*< ELSE >*/
} else {
/*< >*/
dlatdf_(ijob, &zdim, z__, &c__8, rhs, rdsum, rdscal,
ipiv, jpiv);
/*< END IF >*/
}
/* Unpack solution vector(s) */
/*< C( IS, JS ) = RHS( 1 ) >*/
c__[is + js * c_dim1] = rhs[0];
/*< F( IS, JS ) = RHS( 2 ) >*/
f[is + js * f_dim1] = rhs[1];
/* Substitute R(I, J) and L(I, J) into remaining */
/* equation. */
/*< IF( I.GT.1 ) THEN >*/
if (i__ > 1) {
/*< ALPHA = -RHS( 1 ) >*/
alpha = -rhs[0];
/*< >*/
i__2 = is - 1;
daxpy_(&i__2, &alpha, &a[is * a_dim1 + 1], &c__1, &
c__[js * c_dim1 + 1], &c__1);
/*< >*/
i__2 = is - 1;
daxpy_(&i__2, &alpha, &d__[is * d_dim1 + 1], &c__1, &
f[js * f_dim1 + 1], &c__1);
/*< END IF >*/
}
/*< IF( J.LT.Q ) THEN >*/
if (j < q) {
/*< >*/
i__2 = *n - je;
daxpy_(&i__2, &rhs[1], &b[js + (je + 1) * b_dim1],
ldb, &c__[is + (je + 1) * c_dim1], ldc);
/*< >*/
i__2 = *n - je;
daxpy_(&i__2, &rhs[1], &e[js + (je + 1) * e_dim1],
lde, &f[is + (je + 1) * f_dim1], ldf);
/*< END IF >*/
}
/*< ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN >*/
} else if (mb == 1 && nb == 2) {
/* Build a 4-by-4 system Z * x = RHS */
/*< Z( 1, 1 ) = A( IS, IS ) >*/
z__[0] = a[is + is * a_dim1];
/*< Z( 2, 1 ) = ZERO >*/
z__[1] = 0.;
/*< Z( 3, 1 ) = D( IS, IS ) >*/
z__[2] = d__[is + is * d_dim1];
/*< Z( 4, 1 ) = ZERO >*/
z__[3] = 0.;
/*< Z( 1, 2 ) = ZERO >*/
z__[8] = 0.;
/*< Z( 2, 2 ) = A( IS, IS ) >*/
z__[9] = a[is + is * a_dim1];
/*< Z( 3, 2 ) = ZERO >*/
z__[10] = 0.;
/*< Z( 4, 2 ) = D( IS, IS ) >*/
z__[11] = d__[is + is * d_dim1];
/*< Z( 1, 3 ) = -B( JS, JS ) >*/
z__[16] = -b[js + js * b_dim1];
/*< Z( 2, 3 ) = -B( JS, JSP1 ) >*/
z__[17] = -b[js + jsp1 * b_dim1];
/*< Z( 3, 3 ) = -E( JS, JS ) >*/
z__[18] = -e[js + js * e_dim1];
/*< Z( 4, 3 ) = -E( JS, JSP1 ) >*/
z__[19] = -e[js + jsp1 * e_dim1];
/*< Z( 1, 4 ) = -B( JSP1, JS ) >*/
z__[24] = -b[jsp1 + js * b_dim1];
/*< Z( 2, 4 ) = -B( JSP1, JSP1 ) >*/
z__[25] = -b[jsp1 + jsp1 * b_dim1];
/*< Z( 3, 4 ) = ZERO >*/
z__[26] = 0.;
/*< Z( 4, 4 ) = -E( JSP1, JSP1 ) >*/
z__[27] = -e[jsp1 + jsp1 * e_dim1];
/* Set up right hand side(s) */
/*< RHS( 1 ) = C( IS, JS ) >*/
rhs[0] = c__[is + js * c_dim1];
/*< RHS( 2 ) = C( IS, JSP1 ) >*/
rhs[1] = c__[is + jsp1 * c_dim1];
/*< RHS( 3 ) = F( IS, JS ) >*/
rhs[2] = f[is + js * f_dim1];
/*< RHS( 4 ) = F( IS, JSP1 ) >*/
rhs[3] = f[is + jsp1 * f_dim1];
/* Solve Z * x = RHS */
/*< CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) >*/
dgetc2_(&zdim, z__, &c__8, ipiv, jpiv, &ierr);
/*< >*/
if (ierr > 0) {
*info = ierr;
}
/*< IF( IJOB.EQ.0 ) THEN >*/
if (*ijob == 0) {
/*< >*/
dgesc2_(&zdim, z__, &c__8, rhs, ipiv, jpiv, &scaloc);
/*< IF( SCALOC.NE.ONE ) THEN >*/
if (scaloc != 1.) {
/*< DO 60 K = 1, N >*/
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
/*< CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &c__[k * c_dim1 + 1], &
c__1);
/*< CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
/*< 60 CONTINUE >*/
/* L60: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/*< ELSE >*/
} else {
/*< >*/
dlatdf_(ijob, &zdim, z__, &c__8, rhs, rdsum, rdscal,
ipiv, jpiv);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -