📄 zlatrs.c
字号:
if (tjj > smlnum) {
/* abs(A(j,j)) > SMLNUM: */
/*< IF( TJJ.LT.ONE ) THEN >*/
if (tjj < 1.) {
/*< IF( XJ.GT.TJJ*BIGNUM ) THEN >*/
if (xj > tjj * bignum) {
/* Scale x by 1/b(j). */
/*< REC = ONE / XJ >*/
rec = 1. / xj;
/*< CALL ZDSCAL( N, REC, X, 1 ) >*/
zdscal_(n, &rec, &x[1], &c__1);
/*< SCALE = SCALE*REC >*/
*scale *= rec;
/*< XMAX = XMAX*REC >*/
xmax *= rec;
/*< END IF >*/
}
/*< END IF >*/
}
/*< X( J ) = ZLADIV( X( J ), TJJS ) >*/
i__3 = j;
zladiv_(&z__1, &x[j], &tjjs);
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/*< XJ = CABS1( X( J ) ) >*/
i__3 = j;
xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j])
, abs(d__2));
/*< ELSE IF( TJJ.GT.ZERO ) THEN >*/
} else if (tjj > 0.) {
/* 0 < abs(A(j,j)) <= SMLNUM: */
/*< IF( XJ.GT.TJJ*BIGNUM ) THEN >*/
if (xj > tjj * bignum) {
/* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM */
/* to avoid overflow when dividing by A(j,j). */
/*< REC = ( TJJ*BIGNUM ) / XJ >*/
rec = tjj * bignum / xj;
/*< IF( CNORM( J ).GT.ONE ) THEN >*/
if (cnorm[j] > 1.) {
/* Scale by 1/CNORM(j) to avoid overflow when */
/* multiplying x(j) times column j. */
/*< REC = REC / CNORM( J ) >*/
rec /= cnorm[j];
/*< END IF >*/
}
/*< CALL ZDSCAL( N, REC, X, 1 ) >*/
zdscal_(n, &rec, &x[1], &c__1);
/*< SCALE = SCALE*REC >*/
*scale *= rec;
/*< XMAX = XMAX*REC >*/
xmax *= rec;
/*< END IF >*/
}
/*< X( J ) = ZLADIV( X( J ), TJJS ) >*/
i__3 = j;
zladiv_(&z__1, &x[j], &tjjs);
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/*< XJ = CABS1( X( J ) ) >*/
i__3 = j;
xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j])
, abs(d__2));
/*< ELSE >*/
} else {
/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */
/* scale = 0, and compute a solution to A*x = 0. */
/*< DO 100 I = 1, N >*/
i__3 = *n;
for (i__ = 1; i__ <= i__3; ++i__) {
/*< X( I ) = ZERO >*/
i__4 = i__;
x[i__4].r = 0., x[i__4].i = 0.;
/*< 100 CONTINUE >*/
/* L100: */
}
/*< X( J ) = ONE >*/
i__3 = j;
x[i__3].r = 1., x[i__3].i = 0.;
/*< XJ = ONE >*/
xj = 1.;
/*< SCALE = ZERO >*/
*scale = 0.;
/*< XMAX = ZERO >*/
xmax = 0.;
/*< END IF >*/
}
/*< 110 CONTINUE >*/
L110:
/* Scale x if necessary to avoid overflow when adding a */
/* multiple of column j of A. */
/*< IF( XJ.GT.ONE ) THEN >*/
if (xj > 1.) {
/*< REC = ONE / XJ >*/
rec = 1. / xj;
/*< IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN >*/
if (cnorm[j] > (bignum - xmax) * rec) {
/* Scale x by 1/(2*abs(x(j))). */
/*< REC = REC*HALF >*/
rec *= .5;
/*< CALL ZDSCAL( N, REC, X, 1 ) >*/
zdscal_(n, &rec, &x[1], &c__1);
/*< SCALE = SCALE*REC >*/
*scale *= rec;
/*< END IF >*/
}
/*< ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN >*/
} else if (xj * cnorm[j] > bignum - xmax) {
/* Scale x by 1/2. */
/*< CALL ZDSCAL( N, HALF, X, 1 ) >*/
zdscal_(n, &c_b36, &x[1], &c__1);
/*< SCALE = SCALE*HALF >*/
*scale *= .5;
/*< END IF >*/
}
/*< IF( UPPER ) THEN >*/
if (upper) {
/*< IF( J.GT.1 ) THEN >*/
if (j > 1) {
/* Compute the update */
/* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) */
/*< >*/
i__3 = j - 1;
i__4 = j;
z__2.r = -x[i__4].r, z__2.i = -x[i__4].i;
z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
zaxpy_(&i__3, &z__1, &a[j * a_dim1 + 1], &c__1, &x[1],
&c__1);
/*< I = IZAMAX( J-1, X, 1 ) >*/
i__3 = j - 1;
i__ = izamax_(&i__3, &x[1], &c__1);
/*< XMAX = CABS1( X( I ) ) >*/
i__3 = i__;
xmax = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(
&x[i__]), abs(d__2));
/*< END IF >*/
}
/*< ELSE >*/
} else {
/*< IF( J.LT.N ) THEN >*/
if (j < *n) {
/* Compute the update */
/* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) */
/*< >*/
i__3 = *n - j;
i__4 = j;
z__2.r = -x[i__4].r, z__2.i = -x[i__4].i;
z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
zaxpy_(&i__3, &z__1, &a[j + 1 + j * a_dim1], &c__1, &
x[j + 1], &c__1);
/*< I = J + IZAMAX( N-J, X( J+1 ), 1 ) >*/
i__3 = *n - j;
i__ = j + izamax_(&i__3, &x[j + 1], &c__1);
/*< XMAX = CABS1( X( I ) ) >*/
i__3 = i__;
xmax = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(
&x[i__]), abs(d__2));
/*< END IF >*/
}
/*< END IF >*/
}
/*< 120 CONTINUE >*/
/* L120: */
}
/*< ELSE IF( LSAME( TRANS, 'T' ) ) THEN >*/
} else if (lsame_(trans, "T", (ftnlen)1, (ftnlen)1)) {
/* Solve A**T * x = b */
/*< DO 170 J = JFIRST, JLAST, JINC >*/
i__2 = jlast;
i__1 = jinc;
for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
/* Compute x(j) = b(j) - sum A(k,j)*x(k). */
/* k<>j */
/*< XJ = CABS1( X( J ) ) >*/
i__3 = j;
xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j]),
abs(d__2));
/*< USCAL = TSCAL >*/
uscal.r = tscal, uscal.i = 0.;
/*< REC = ONE / MAX( XMAX, ONE ) >*/
rec = 1. / max(xmax,1.);
/*< IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN >*/
if (cnorm[j] > (bignum - xj) * rec) {
/* If x(j) could overflow, scale x by 1/(2*XMAX). */
/*< REC = REC*HALF >*/
rec *= .5;
/*< IF( NOUNIT ) THEN >*/
if (nounit) {
/*< TJJS = A( J, J )*TSCAL >*/
i__3 = j + j * a_dim1;
z__1.r = tscal * a[i__3].r, z__1.i = tscal * a[i__3]
.i;
tjjs.r = z__1.r, tjjs.i = z__1.i;
/*< ELSE >*/
} else {
/*< TJJS = TSCAL >*/
tjjs.r = tscal, tjjs.i = 0.;
/*< END IF >*/
}
/*< TJJ = CABS1( TJJS ) >*/
tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs),
abs(d__2));
/*< IF( TJJ.GT.ONE ) THEN >*/
if (tjj > 1.) {
/* Divide by A(j,j) when scaling x if A(j,j) > 1. */
/*< REC = MIN( ONE, REC*TJJ ) >*/
/* Computing MIN */
d__1 = 1., d__2 = rec * tjj;
rec = min(d__1,d__2);
/*< USCAL = ZLADIV( USCAL, TJJS ) >*/
zladiv_(&z__1, &uscal, &tjjs);
uscal.r = z__1.r, uscal.i = z__1.i;
/*< END IF >*/
}
/*< IF( REC.LT.ONE ) THEN >*/
if (rec < 1.) {
/*< CALL ZDSCAL( N, REC, X, 1 ) >*/
zdscal_(n, &rec, &x[1], &c__1);
/*< SCALE = SCALE*REC >*/
*scale *= rec;
/*< XMAX = XMAX*REC >*/
xmax *= rec;
/*< END IF >*/
}
/*< END IF >*/
}
/*< CSUMJ = ZERO >*/
csumj.r = 0., csumj.i = 0.;
/*< IF( USCAL.EQ.DCMPLX( ONE ) ) THEN >*/
if (uscal.r == 1. && uscal.i == 0.) {
/* If the scaling needed for A in the dot product is 1, */
/* call ZDOTU to perform the dot product. */
/*< IF( UPPER ) THEN >*/
if (upper) {
/*< CSUMJ = ZDOTU( J-1, A( 1, J ), 1, X, 1 ) >*/
i__3 = j - 1;
zdotu_(&z__1, &i__3, &a[j * a_dim1 + 1], &c__1, &x[1],
&c__1);
csumj.r = z__1.r, csumj.i = z__1.i;
/*< ELSE IF( J.LT.N ) THEN >*/
} else if (j < *n) {
/*< CSUMJ = ZDOTU( N-J, A( J+1, J ), 1, X( J+1 ), 1 ) >*/
i__3 = *n - j;
zdotu_(&z__1, &i__3, &a[j + 1 + j * a_dim1], &c__1, &
x[j + 1], &c__1);
csumj.r = z__1.r, csumj.i = z__1.i;
/*< END IF >*/
}
/*< ELSE >*/
} else {
/* Otherwise, use in-line code for the dot product. */
/*< IF( UPPER ) THEN >*/
if (upper) {
/*< DO 130 I = 1, J - 1 >*/
i__3 = j - 1;
for (i__ = 1; i__ <= i__3; ++i__) {
/*< CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I ) >*/
i__4 = i__ + j * a_dim1;
z__3.r = a[i__4].r * uscal.r - a[i__4].i *
uscal.i, z__3.i = a[i__4].r * uscal.i + a[
i__4].i * uscal.r;
i__5 = i__;
z__2.r = z__3.r * x[i__5].r - z__3.i * x[i__5].i,
z__2.i = z__3.r * x[i__5].i + z__3.i * x[
i__5].r;
z__1.r = csumj.r + z__2.r, z__1.i = csumj.i +
z__2.i;
csumj.r = z__1.r, csumj.i = z__1.i;
/*< 130 CONTINUE >*/
/* L130: */
}
/*< ELSE IF( J.LT.N ) THEN >*/
} else if (j < *n) {
/*< DO 140 I = J + 1, N >*/
i__3 = *n;
for (i__ = j + 1; i__ <= i__3; ++i__) {
/*< CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I ) >*/
i__4 = i__ + j * a_dim1;
z__3.r = a[i__4].r * uscal.r - a[i__4].i *
uscal.i, z__3.i = a[i__4].r * uscal.i + a[
i__4].i * uscal.r;
i__5 = i__;
z__2.r = z__3.r * x[i__5].r - z__3.i * x[i__5].i,
z__2.i = z__3.r * x[i__5].i + z__3.i * x[
i__5].r;
z__1.r = csumj.r + z__2.r, z__1.i = csumj.i +
z__2.i;
csumj.r = z__1.r, csumj.i = z__1.i;
/*< 140 CONTINUE >*/
/* L140: */
}
/*< END IF >*/
}
/*< END IF >*/
}
/*< IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN >*/
z__1.r = tscal, z__1.i = 0.;
if (uscal.r == z__1.r && uscal.i == z__1.i) {
/* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) */
/* was not used to scale the dotproduct. */
/*< X( J ) = X( J ) - CSUMJ >*/
i__3 = j;
i__4 = j;
z__1.r = x[i__4].r - csumj.r, z__1.i = x[i__4].i -
csumj.i;
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/*< XJ = CABS1( X( J ) ) >*/
i__3 = j;
xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j])
, abs(d__2));
/*< IF( NOUNIT ) THEN >*/
if (nounit) {
/*< TJJS = A( J, J )*TSCAL >*/
i__3 = j + j * a_dim1;
z__1.r = tscal * a[i__3].r, z__1.i = tscal * a[i__3]
.i;
tjjs.r = z__1.r, tjjs.i = z__1.i;
/*< ELSE >*/
} else {
/*< TJJS = TSCAL >*/
tjjs.r = tscal, tjjs.i = 0.;
/*< >*/
if (tscal == 1.) {
goto L160;
}
/*< END IF >*/
}
/* Compute x(j) = x(j) / A(j,j), scaling if necessary. */
/*< TJJ = CABS1( TJJS ) >*/
tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs),
abs(d__2));
/*< IF( TJJ.GT.SMLNUM ) THEN >*/
if (tjj > smlnum) {
/* abs(A(j,j)) > SMLNUM: */
/*< IF( TJJ.LT.ONE ) THEN >*/
if (tjj < 1.) {
/*< IF( XJ.GT.TJJ*BIGNUM ) THEN >*/
if (xj > tjj * bignum) {
/* Scale X by 1/abs(x(j)). */
/*< REC = ONE / XJ >*/
rec = 1. / xj;
/*< CALL ZDSCAL( N, REC, X, 1 ) >*/
zdscal_(n, &rec, &x[1], &c__1);
/*< SCALE = SCALE*REC >*/
*scale *= rec;
/*< XMAX = XMAX*REC >*/
xmax *= rec;
/*< END IF >*/
}
/*< END IF >*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -