📄 zlatrs.c
字号:
}
/*< 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;
/*< 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. */
/*< REC = ( TJJ*BIGNUM ) / XJ >*/
rec = tjj * bignum / 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 >*/
}
/*< 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;
/*< ELSE >*/
} else {
/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */
/* scale = 0 and compute a solution to A**T *x = 0. */
/*< DO 150 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.;
/*< 150 CONTINUE >*/
/* L150: */
}
/*< X( J ) = ONE >*/
i__3 = j;
x[i__3].r = 1., x[i__3].i = 0.;
/*< SCALE = ZERO >*/
*scale = 0.;
/*< XMAX = ZERO >*/
xmax = 0.;
/*< END IF >*/
}
/*< 160 CONTINUE >*/
L160:
/*< ELSE >*/
;
} else {
/* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot */
/* product has already been divided by 1/A(j,j). */
/*< X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ >*/
i__3 = j;
zladiv_(&z__2, &x[j], &tjjs);
z__1.r = z__2.r - csumj.r, z__1.i = z__2.i - csumj.i;
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/*< END IF >*/
}
/*< XMAX = MAX( XMAX, CABS1( X( J ) ) ) >*/
/* Computing MAX */
i__3 = j;
d__3 = xmax, d__4 = (d__1 = x[i__3].r, abs(d__1)) + (d__2 =
d_imag(&x[j]), abs(d__2));
xmax = max(d__3,d__4);
/*< 170 CONTINUE >*/
/* L170: */
}
/*< ELSE >*/
} else {
/* Solve A**H * x = b */
/*< DO 220 J = JFIRST, JLAST, JINC >*/
i__1 = jlast;
i__2 = jinc;
for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
/* 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 = DCONJG( A( J, J ) )*TSCAL >*/
d_cnjg(&z__2, &a[j + j * a_dim1]);
z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.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 ZDOTC to perform the dot product. */
/*< IF( UPPER ) THEN >*/
if (upper) {
/*< CSUMJ = ZDOTC( J-1, A( 1, J ), 1, X, 1 ) >*/
i__3 = j - 1;
zdotc_(&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 = ZDOTC( N-J, A( J+1, J ), 1, X( J+1 ), 1 ) >*/
i__3 = *n - j;
zdotc_(&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 180 I = 1, J - 1 >*/
i__3 = j - 1;
for (i__ = 1; i__ <= i__3; ++i__) {
/*< >*/
d_cnjg(&z__4, &a[i__ + j * a_dim1]);
z__3.r = z__4.r * uscal.r - z__4.i * uscal.i,
z__3.i = z__4.r * uscal.i + z__4.i *
uscal.r;
i__4 = i__;
z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i,
z__2.i = z__3.r * x[i__4].i + z__3.i * x[
i__4].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;
/*< 180 CONTINUE >*/
/* L180: */
}
/*< ELSE IF( J.LT.N ) THEN >*/
} else if (j < *n) {
/*< DO 190 I = J + 1, N >*/
i__3 = *n;
for (i__ = j + 1; i__ <= i__3; ++i__) {
/*< >*/
d_cnjg(&z__4, &a[i__ + j * a_dim1]);
z__3.r = z__4.r * uscal.r - z__4.i * uscal.i,
z__3.i = z__4.r * uscal.i + z__4.i *
uscal.r;
i__4 = i__;
z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i,
z__2.i = z__3.r * x[i__4].i + z__3.i * x[
i__4].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;
/*< 190 CONTINUE >*/
/* L190: */
}
/*< 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 = DCONJG( A( J, J ) )*TSCAL >*/
d_cnjg(&z__2, &a[j + j * a_dim1]);
z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.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 L210;
}
/*< 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 >*/
}
/*< 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;
/*< 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. */
/*< REC = ( TJJ*BIGNUM ) / XJ >*/
rec = tjj * bignum / 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 >*/
}
/*< 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;
/*< ELSE >*/
} else {
/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */
/* scale = 0 and compute a solution to A**H *x = 0. */
/*< DO 200 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.;
/*< 200 CONTINUE >*/
/* L200: */
}
/*< X( J ) = ONE >*/
i__3 = j;
x[i__3].r = 1., x[i__3].i = 0.;
/*< SCALE = ZERO >*/
*scale = 0.;
/*< XMAX = ZERO >*/
xmax = 0.;
/*< END IF >*/
}
/*< 210 CONTINUE >*/
L210:
/*< ELSE >*/
;
} else {
/* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot */
/* product has already been divided by 1/A(j,j). */
/*< X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ >*/
i__3 = j;
zladiv_(&z__2, &x[j], &tjjs);
z__1.r = z__2.r - csumj.r, z__1.i = z__2.i - csumj.i;
x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/*< END IF >*/
}
/*< XMAX = MAX( XMAX, CABS1( X( J ) ) ) >*/
/* Computing MAX */
i__3 = j;
d__3 = xmax, d__4 = (d__1 = x[i__3].r, abs(d__1)) + (d__2 =
d_imag(&x[j]), abs(d__2));
xmax = max(d__3,d__4);
/*< 220 CONTINUE >*/
/* L220: */
}
/*< END IF >*/
}
/*< SCALE = SCALE / TSCAL >*/
*scale /= tscal;
/*< END IF >*/
}
/* Scale the column norms by 1/TSCAL for return. */
/*< IF( TSCAL.NE.ONE ) THEN >*/
if (tscal != 1.) {
/*< CALL DSCAL( N, ONE / TSCAL, CNORM, 1 ) >*/
d__1 = 1. / tscal;
dscal_(n, &d__1, &cnorm[1], &c__1);
/*< END IF >*/
}
/*< RETURN >*/
return 0;
/* End of ZLATRS */
/*< END >*/
} /* zlatrs_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -