zlatrs.c
来自「算断裂的」· C语言 代码 · 共 1,187 行 · 第 1/3 页
C
1,187 行
d__1 = xbnd, d__2 = min(1.,tjj) * grow;
xbnd = min(d__1,d__2);
} else {
/* M(j) could overflow, set XBND to 0. */
xbnd = 0.;
}
if (tjj + CNORM(j) >= smlnum) {
/* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,
j)) ) */
grow *= tjj / (tjj + CNORM(j));
} else {
/* G(j) could overflow, set GROW to 0. */
grow = 0.;
}
/* L40: */
}
grow = xbnd;
} else {
/* A is unit triangular.
Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...
,n}.
Computing MIN */
d__1 = 1., d__2 = .5 / max(xbnd,smlnum);
grow = min(d__1,d__2);
i__2 = jlast;
i__1 = jinc;
for (j = jfirst; jinc < 0 ? j >= jlast : j <= jlast; j += jinc) {
/* Exit the loop if the growth factor is too smal
l. */
if (grow <= smlnum) {
goto L60;
}
/* G(j) = G(j-1)*( 1 + CNORM(j) ) */
grow *= 1. / (CNORM(j) + 1.);
/* L50: */
}
}
L60:
;
} else {
/* Compute the growth in A**T * x = b or A**H * x = b. */
if (upper) {
jfirst = 1;
jlast = *n;
jinc = 1;
} else {
jfirst = *n;
jlast = 1;
jinc = -1;
}
if (tscal != 1.) {
grow = 0.;
goto L90;
}
if (nounit) {
/* A is non-unit triangular.
Compute GROW = 1/G(j) and XBND = 1/M(j).
Initially, M(0) = max{x(i), i=1,...,n}. */
grow = .5 / max(xbnd,smlnum);
xbnd = grow;
i__1 = jlast;
i__2 = jinc;
for (j = jfirst; jinc < 0 ? j >= jlast : j <= jlast; j += jinc) {
/* Exit the loop if the growth factor is too smal
l. */
if (grow <= smlnum) {
goto L90;
}
/* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
*/
xj = CNORM(j) + 1.;
/* Computing MIN */
d__1 = grow, d__2 = xbnd / xj;
grow = min(d__1,d__2);
i__3 = j + j * a_dim1;
tjjs.r = A(j,j).r, tjjs.i = A(j,j).i;
tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(
d__2));
if (tjj >= smlnum) {
/* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(
j,j)) */
if (xj > tjj) {
xbnd *= tjj / xj;
}
} else {
/* M(j) could overflow, set XBND to 0. */
xbnd = 0.;
}
/* L70: */
}
grow = min(grow,xbnd);
} else {
/* A is unit triangular.
Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...
,n}.
Computing MIN */
d__1 = 1., d__2 = .5 / max(xbnd,smlnum);
grow = min(d__1,d__2);
i__2 = jlast;
i__1 = jinc;
for (j = jfirst; jinc < 0 ? j >= jlast : j <= jlast; j += jinc) {
/* Exit the loop if the growth factor is too smal
l. */
if (grow <= smlnum) {
goto L90;
}
/* G(j) = ( 1 + CNORM(j) )*G(j-1) */
xj = CNORM(j) + 1.;
grow /= xj;
/* L80: */
}
}
L90:
;
}
if (grow * tscal > smlnum) {
/* Use the Level 2 BLAS solve if the reciprocal of the bound on
elements of X is not too small. */
ztrsv_(uplo, trans, diag, n, &A(1,1), lda, &X(1), &c__1);
} else {
/* Use a Level 1 BLAS solve, scaling intermediate results. */
if (xmax > bignum * .5) {
/* Scale X so that its components are less than or equal
to
BIGNUM in absolute value. */
*scale = bignum * .5 / xmax;
zdscal_(n, scale, &X(1), &c__1);
xmax = bignum;
} else {
xmax *= 2.;
}
if (notran) {
/* Solve A * x = b */
i__1 = jlast;
i__2 = jinc;
for (j = jfirst; jinc < 0 ? j >= jlast : j <= jlast; j += jinc) {
/* Compute x(j) = b(j) / A(j,j), scaling x if nec
essary. */
i__3 = j;
xj = (d__1 = X(j).r, abs(d__1)) + (d__2 = d_imag(&X(j)),
abs(d__2));
if (nounit) {
i__3 = j + j * a_dim1;
z__1.r = tscal * A(j,j).r, z__1.i = tscal * A(j,j).i;
tjjs.r = z__1.r, tjjs.i = z__1.i;
} else {
tjjs.r = tscal, tjjs.i = 0.;
if (tscal == 1.) {
goto L110;
}
}
tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(
d__2));
if (tjj > smlnum) {
/* abs(A(j,j)) > SMLNUM: */
if (tjj < 1.) {
if (xj > tjj * bignum) {
/* Scale x by 1/b(j). */
rec = 1. / xj;
zdscal_(n, &rec, &X(1), &c__1);
*scale *= rec;
xmax *= rec;
}
}
i__3 = j;
zladiv_(&z__1, &X(j), &tjjs);
X(j).r = z__1.r, X(j).i = z__1.i;
i__3 = j;
xj = (d__1 = X(j).r, abs(d__1)) + (d__2 = d_imag(&X(j))
, abs(d__2));
} else if (tjj > 0.) {
/* 0 < abs(A(j,j)) <= SMLNUM: */
if (xj > tjj * bignum) {
/* Scale x by (1/abs(x(j)))*abs(
A(j,j))*BIGNUM
to avoid overflow when dividi
ng by A(j,j). */
rec = tjj * bignum / xj;
if (CNORM(j) > 1.) {
/* Scale by 1/CNORM(j) to
avoid overflow when
multiplying x(j) times
column j. */
rec /= CNORM(j);
}
zdscal_(n, &rec, &X(1), &c__1);
*scale *= rec;
xmax *= rec;
}
i__3 = j;
zladiv_(&z__1, &X(j), &tjjs);
X(j).r = z__1.r, X(j).i = z__1.i;
i__3 = j;
xj = (d__1 = X(j).r, abs(d__1)) + (d__2 = d_imag(&X(j))
, abs(d__2));
} else {
/* A(j,j) = 0: Set x(1:n) = 0, x(j) =
1, and
scale = 0, and compute a solution to
A*x = 0. */
i__3 = *n;
for (i = 1; i <= *n; ++i) {
i__4 = i;
X(i).r = 0., X(i).i = 0.;
/* L100: */
}
i__3 = j;
X(j).r = 1., X(j).i = 0.;
xj = 1.;
*scale = 0.;
xmax = 0.;
}
L110:
/* Scale x if necessary to avoid overflow when ad
ding a
multiple of column j of A. */
if (xj > 1.) {
rec = 1. / xj;
if (CNORM(j) > (bignum - xmax) * rec) {
/* Scale x by 1/(2*abs(x(j))). */
rec *= .5;
zdscal_(n, &rec, &X(1), &c__1);
*scale *= rec;
}
} else if (xj * CNORM(j) > bignum - xmax) {
/* Scale x by 1/2. */
zdscal_(n, &c_b36, &X(1), &c__1);
*scale *= .5;
}
if (upper) {
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(j).r, z__2.i = -X(j).i;
z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
zaxpy_(&i__3, &z__1, &A(1,j), &c__1, &X(1),
&c__1);
i__3 = j - 1;
i = izamax_(&i__3, &X(1), &c__1);
i__3 = i;
xmax = (d__1 = X(i).r, abs(d__1)) + (d__2 = d_imag(
&X(i)), abs(d__2));
}
} else {
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(j).r, z__2.i = -X(j).i;
z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
zaxpy_(&i__3, &z__1, &A(j+1,j), &c__1, &
X(j + 1), &c__1);
i__3 = *n - j;
i = j + izamax_(&i__3, &X(j + 1), &c__1);
i__3 = i;
xmax = (d__1 = X(i).r, abs(d__1)) + (d__2 = d_imag(
&X(i)), abs(d__2));
}
}
/* L120: */
}
} else if (lsame_(trans, "T")) {
/* Solve A**T * x = b */
i__2 = jlast;
i__1 = jinc;
for (j = jfirst; jinc < 0 ? j >= jlast : j <= jlast; j += jinc) {
/* Compute x(j) = b(j) - sum A(k,j)*x(k).
k<>j */
i__3 = j;
xj = (d__1 = X(j).r, abs(d__1)) + (d__2 = d_imag(&X(j)),
abs(d__2));
uscal.r = tscal, uscal.i = 0.;
rec = 1. / max(xmax,1.);
if (CNORM(j) > (bignum - xj) * rec) {
/* If x(j) could overflow, scale x by 1/(2
*XMAX). */
rec *= .5;
if (nounit) {
i__3 = j + j * a_dim1;
z__1.r = tscal * A(j,j).r, z__1.i = tscal * A(j,j)
.i;
tjjs.r = z__1.r, tjjs.i = z__1.i;
} else {
tjjs.r = tscal, tjjs.i = 0.;
}
tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs),
abs(d__2));
if (tjj > 1.) {
/* Divide by A(j,j) when scaling
x if A(j,j) > 1.
Computing MIN */
d__1 = 1., d__2 = rec * tjj;
rec = min(d__1,d__2);
zladiv_(&z__1, &uscal, &tjjs);
uscal.r = z__1.r, uscal.i = z__1.i;
}
if (rec < 1.) {
zdscal_(n, &rec, &X(1), &c__1);
*scale *= rec;
xmax *= rec;
}
}
csumj.r = 0., csumj.i = 0.;
if (uscal.r == 1. && uscal.i == 0.) {
/* If the scaling needed for A in the dot
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?