📄 dlatrs.c
字号:
xmax *= rec;
/*< END IF >*/
}
/*< X( J ) = X( J ) / TJJS >*/
x[j] /= tjjs;
/*< XJ = ABS( X( J ) ) >*/
xj = (d__1 = x[j], abs(d__1));
/*< 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 90 I = 1, N >*/
i__3 = *n;
for (i__ = 1; i__ <= i__3; ++i__) {
/*< X( I ) = ZERO >*/
x[i__] = 0.;
/*< 90 CONTINUE >*/
/* L90: */
}
/*< X( J ) = ONE >*/
x[j] = 1.;
/*< XJ = ONE >*/
xj = 1.;
/*< SCALE = ZERO >*/
*scale = 0.;
/*< XMAX = ZERO >*/
xmax = 0.;
/*< END IF >*/
}
/*< 100 CONTINUE >*/
L100:
/* 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 DSCAL( N, REC, X, 1 ) >*/
dscal_(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 DSCAL( N, HALF, X, 1 ) >*/
dscal_(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;
d__1 = -x[j] * tscal;
daxpy_(&i__3, &d__1, &a[j * a_dim1 + 1], &c__1, &x[1],
&c__1);
/*< I = IDAMAX( J-1, X, 1 ) >*/
i__3 = j - 1;
i__ = idamax_(&i__3, &x[1], &c__1);
/*< XMAX = ABS( X( I ) ) >*/
xmax = (d__1 = x[i__], abs(d__1));
/*< 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;
d__1 = -x[j] * tscal;
daxpy_(&i__3, &d__1, &a[j + 1 + j * a_dim1], &c__1, &
x[j + 1], &c__1);
/*< I = J + IDAMAX( N-J, X( J+1 ), 1 ) >*/
i__3 = *n - j;
i__ = j + idamax_(&i__3, &x[j + 1], &c__1);
/*< XMAX = ABS( X( I ) ) >*/
xmax = (d__1 = x[i__], abs(d__1));
/*< END IF >*/
}
/*< END IF >*/
}
/*< 110 CONTINUE >*/
/* L110: */
}
/*< ELSE >*/
} else {
/* Solve A' * x = b */
/*< DO 160 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 = ABS( X( J ) ) >*/
xj = (d__1 = x[j], abs(d__1));
/*< USCAL = TSCAL >*/
uscal = tscal;
/*< 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 >*/
tjjs = a[j + j * a_dim1] * tscal;
/*< ELSE >*/
} else {
/*< TJJS = TSCAL >*/
tjjs = tscal;
/*< END IF >*/
}
/*< TJJ = ABS( TJJS ) >*/
tjj = abs(tjjs);
/*< 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 = USCAL / TJJS >*/
uscal /= tjjs;
/*< END IF >*/
}
/*< IF( REC.LT.ONE ) THEN >*/
if (rec < 1.) {
/*< CALL DSCAL( N, REC, X, 1 ) >*/
dscal_(n, &rec, &x[1], &c__1);
/*< SCALE = SCALE*REC >*/
*scale *= rec;
/*< XMAX = XMAX*REC >*/
xmax *= rec;
/*< END IF >*/
}
/*< END IF >*/
}
/*< SUMJ = ZERO >*/
sumj = 0.;
/*< IF( USCAL.EQ.ONE ) THEN >*/
if (uscal == 1.) {
/* If the scaling needed for A in the dot product is 1, */
/* call DDOT to perform the dot product. */
/*< IF( UPPER ) THEN >*/
if (upper) {
/*< SUMJ = DDOT( J-1, A( 1, J ), 1, X, 1 ) >*/
i__3 = j - 1;
sumj = ddot_(&i__3, &a[j * a_dim1 + 1], &c__1, &x[1],
&c__1);
/*< ELSE IF( J.LT.N ) THEN >*/
} else if (j < *n) {
/*< SUMJ = DDOT( N-J, A( J+1, J ), 1, X( J+1 ), 1 ) >*/
i__3 = *n - j;
sumj = ddot_(&i__3, &a[j + 1 + j * a_dim1], &c__1, &x[
j + 1], &c__1);
/*< END IF >*/
}
/*< ELSE >*/
} else {
/* Otherwise, use in-line code for the dot product. */
/*< IF( UPPER ) THEN >*/
if (upper) {
/*< DO 120 I = 1, J - 1 >*/
i__3 = j - 1;
for (i__ = 1; i__ <= i__3; ++i__) {
/*< SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I ) >*/
sumj += a[i__ + j * a_dim1] * uscal * x[i__];
/*< 120 CONTINUE >*/
/* L120: */
}
/*< ELSE IF( J.LT.N ) THEN >*/
} else if (j < *n) {
/*< DO 130 I = J + 1, N >*/
i__3 = *n;
for (i__ = j + 1; i__ <= i__3; ++i__) {
/*< SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I ) >*/
sumj += a[i__ + j * a_dim1] * uscal * x[i__];
/*< 130 CONTINUE >*/
/* L130: */
}
/*< END IF >*/
}
/*< END IF >*/
}
/*< IF( USCAL.EQ.TSCAL ) THEN >*/
if (uscal == tscal) {
/* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j) */
/* was not used to scale the dotproduct. */
/*< X( J ) = X( J ) - SUMJ >*/
x[j] -= sumj;
/*< XJ = ABS( X( J ) ) >*/
xj = (d__1 = x[j], abs(d__1));
/*< IF( NOUNIT ) THEN >*/
if (nounit) {
/*< TJJS = A( J, J )*TSCAL >*/
tjjs = a[j + j * a_dim1] * tscal;
/*< ELSE >*/
} else {
/*< TJJS = TSCAL >*/
tjjs = tscal;
/*< >*/
if (tscal == 1.) {
goto L150;
}
/*< END IF >*/
}
/* Compute x(j) = x(j) / A(j,j), scaling if necessary. */
/*< TJJ = ABS( TJJS ) >*/
tjj = abs(tjjs);
/*< 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 DSCAL( N, REC, X, 1 ) >*/
dscal_(n, &rec, &x[1], &c__1);
/*< SCALE = SCALE*REC >*/
*scale *= rec;
/*< XMAX = XMAX*REC >*/
xmax *= rec;
/*< END IF >*/
}
/*< END IF >*/
}
/*< X( J ) = X( J ) / TJJS >*/
x[j] /= tjjs;
/*< 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 DSCAL( N, REC, X, 1 ) >*/
dscal_(n, &rec, &x[1], &c__1);
/*< SCALE = SCALE*REC >*/
*scale *= rec;
/*< XMAX = XMAX*REC >*/
xmax *= rec;
/*< END IF >*/
}
/*< X( J ) = X( J ) / TJJS >*/
x[j] /= tjjs;
/*< 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 140 I = 1, N >*/
i__3 = *n;
for (i__ = 1; i__ <= i__3; ++i__) {
/*< X( I ) = ZERO >*/
x[i__] = 0.;
/*< 140 CONTINUE >*/
/* L140: */
}
/*< X( J ) = ONE >*/
x[j] = 1.;
/*< SCALE = ZERO >*/
*scale = 0.;
/*< XMAX = ZERO >*/
xmax = 0.;
/*< END IF >*/
}
/*< 150 CONTINUE >*/
L150:
/*< ELSE >*/
;
} else {
/* Compute x(j) := x(j) / A(j,j) - sumj if the dot */
/* product has already been divided by 1/A(j,j). */
/*< X( J ) = X( J ) / TJJS - SUMJ >*/
x[j] = x[j] / tjjs - sumj;
/*< END IF >*/
}
/*< XMAX = MAX( XMAX, ABS( X( J ) ) ) >*/
/* Computing MAX */
d__2 = xmax, d__3 = (d__1 = x[j], abs(d__1));
xmax = max(d__2,d__3);
/*< 160 CONTINUE >*/
/* L160: */
}
/*< 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 DLATRS */
/*< END >*/
} /* dlatrs_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -