zlatrs.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 913 行 · 第 1/3 页
C
913 行
/* Compute x(j) = b(j) - sum A(k,j)*x(k). */
/* k<>j */
xj = abs(x[j].r) + abs(x[j].i);
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__1 = j + j * *lda;
tjjs.r = tscal * a[i__1].r, tjjs.i = tscal * a[i__1].i;
} else {
tjjs.r = tscal, tjjs.i = 0.;
}
tjj = abs(tjjs.r) + abs(tjjs.i);
if (tjj > 1.) {
/* Divide by A(j,j) when scaling x if A(j,j) > 1. */
rec = min(1., rec * tjj);
zladiv_(&uscal, &uscal, &tjjs);
}
if (rec < 1.) {
zdscal_(n, &rec, x, &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 product is 1, */
/* call ZDOTU to perform the dot product. */
if (upper) {
zdotu_(&csumj, &j, &a[j * *lda], &c__1, x, &c__1);
} else if (j < *n - 1) {
i__1 = *n - j - 1;
zdotu_(&csumj, &i__1, &a[j + 1 + j * *lda], &c__1, &x[j + 1], &c__1);
}
} else {
/* Otherwise, use in-line code for the dot product. */
if (upper) {
for (i = 0; i < j - 1; ++i) {
i__1 = i + j * *lda;
z__1.r = a[i__1].r * uscal.r - a[i__1].i * uscal.i,
z__1.i = a[i__1].r * uscal.i + a[i__1].i * uscal.r;
csumj.r += z__1.r * x[i].r - z__1.i * x[i].i,
csumj.i += z__1.r * x[i].i + z__1.i * x[i].r;
}
} else if (j < *n - 1) {
for (i = j + 1; i < *n; ++i) {
i__1 = i + j * *lda;
z__1.r = a[i__1].r * uscal.r - a[i__1].i * uscal.i,
z__1.i = a[i__1].r * uscal.i + a[i__1].i * uscal.r;
csumj.r += z__1.r * x[i].r - z__1.i * x[i].i,
csumj.i += z__1.r * x[i].i + z__1.i * x[i].r;
}
}
}
if (uscal.r == tscal && uscal.i == 0.) {
/* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) */
/* was not used to scale the dotproduct. */
x[j].r = x[j].r - csumj.r, x[j].i = x[j].i - csumj.i;
xj = abs(x[j].r) + abs(x[j].i);
if (nounit) {
i__1 = j + j * *lda;
tjjs.r = tscal * a[i__1].r, tjjs.i = tscal * a[i__1].i;
} else {
tjjs.r = tscal, tjjs.i = 0.;
if (tscal == 1.) {
goto L160;
}
}
/* Compute x(j) = x(j) / A(j,j), scaling if necessary. */
tjj = abs(tjjs.r) + abs(tjjs.i);
if (tjj > smlnum) {
/* abs(A(j,j)) > SMLNUM: */
if (tjj < 1.) {
if (xj > tjj * bignum) {
/* Scale X by 1/abs(x(j)). */
rec = 1. / xj;
zdscal_(n, &rec, x, &c__1);
*scale *= rec;
xmax *= rec;
}
}
zladiv_(&x[j], &x[j], &tjjs);
} 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. */
rec = tjj * bignum / xj;
zdscal_(n, &rec, x, &c__1);
*scale *= rec;
xmax *= rec;
}
zladiv_(&x[j], &x[j], &tjjs);
} 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. */
for (i = 0; i < *n; ++i) {
x[i].r = 0., x[i].i = 0.;
}
x[j].r = 1., x[j].i = 0.;
*scale = 0.;
xmax = 0.;
}
L160:
;
} else {
/* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot */
/* product has already been divided by 1/A(j,j). */
zladiv_(&x[j], &x[j], &tjjs);
}
xmax = max(xmax, abs(x[j].r) + abs(x[j].i));
}
} else {
/* Solve A**H * x = b */
for (j = jfirst; j != jlast; j += jinc) {
/* Compute x(j) = b(j) - sum A(k,j)*x(k). */
/* k<>j */
xj = abs(x[j].r) + abs(x[j].i);
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__1 = j + j * *lda;
tjjs.r = tscal * a[i__1].r,
tjjs.i = - tscal * a[i__1].i;
} else {
tjjs.r = tscal, tjjs.i = 0.;
}
tjj = abs(tjjs.r) + abs(tjjs.i);
if (tjj > 1.) {
/* Divide by A(j,j) when scaling x if A(j,j) > 1. */
rec = min(1., rec * tjj);
zladiv_(&uscal, &uscal, &tjjs);
}
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 product is 1, */
/* call ZDOTC to perform the dot product. */
if (upper) {
zdotc_(&csumj, &j, &a[j**lda], &c__1, x, &c__1);
} else if (j < *n - 1) {
i__1 = *n - j - 1;
zdotc_(&csumj, &i__1, &a[j+1+j**lda], &c__1, &x[j+1], &c__1);
}
} else {
/* Otherwise, use in-line code for the dot product. */
if (upper) {
for (i = 0; i < j - 1; ++i) {
i__1 = i + j * *lda;
z__1.r = a[i__1].r * uscal.r + a[i__1].i * uscal.i,
z__1.i = a[i__1].r * uscal.i - a[i__1].i * uscal.r;
csumj.r += z__1.r * x[i].r - z__1.i * x[i].i,
csumj.i += z__1.r * x[i].i + z__1.i * x[i].r;
}
} else if (j < *n - 1) {
for (i = j + 1; i < *n; ++i) {
i__1 = i + j * *lda;
z__1.r = a[i__1].r * uscal.r + a[i__1].i * uscal.i,
z__1.i = a[i__1].r * uscal.i - a[i__1].i * uscal.r;
csumj.r += z__1.r * x[i].r - z__1.i * x[i].i,
csumj.i += z__1.r * x[i].i + z__1.i * x[i].r;
}
}
}
if (uscal.r == tscal && uscal.i == 0.) {
/* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) */
/* was not used to scale the dotproduct. */
x[j].r = x[j].r - csumj.r, x[j].i = x[j].i - csumj.i;
xj = abs(x[j].r) + abs(x[j].i);
if (nounit) {
i__1 = j + j * *lda;
tjjs.r = tscal * a[i__1].r,
tjjs.i = - tscal * a[i__1].i;
} else {
tjjs.r = tscal, tjjs.i = 0.;
if (tscal == 1.) {
goto L210;
}
}
/* Compute x(j) = x(j) / A(j,j), scaling if necessary. */
tjj = abs(tjjs.r) + abs(tjjs.i);
if (tjj > smlnum) {
/* abs(A(j,j)) > SMLNUM: */
if (tjj < 1.) {
if (xj > tjj * bignum) {
/* Scale X by 1/abs(x(j)). */
rec = 1. / xj;
zdscal_(n, &rec, x, &c__1);
*scale *= rec;
xmax *= rec;
}
}
zladiv_(&x[j], &x[j], &tjjs);
} 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. */
rec = tjj * bignum / xj;
zdscal_(n, &rec, x, &c__1);
*scale *= rec;
xmax *= rec;
}
zladiv_(&x[j], &x[j], &tjjs);
} 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. */
for (i = 0; i < *n; ++i) {
x[i].r = 0., x[i].i = 0.;
}
x[j].r = 1., x[j].i = 0.;
*scale = 0.;
xmax = 0.;
}
L210:
;
} else {
/* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot */
/* product has already been divided by 1/A(j,j). */
zladiv_(&x[j], &x[j], &tjjs);
x[j].r -= csumj.r, x[j].i -= csumj.i;
}
xmax = max(xmax, abs(x[j].r) + abs(x[j].i));
}
}
*scale /= tscal;
}
/* Scale the column norms by 1/TSCAL for return. */
if (tscal != 1.) {
d__1 = 1. / tscal;
dscal_(n, &d__1, cnorm, &c__1);
}
} /* zlatrs_ */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?