dlatrs.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*<       IF( TMAX.LE.BIGNUM ) THEN >*/
    if (tmax <= bignum) {
/*<          TSCAL = ONE >*/
        tscal = 1.;
/*<       ELSE >*/
    } else {
/*<          TSCAL = ONE / ( SMLNUM*TMAX ) >*/
        tscal = 1. / (smlnum * tmax);
/*<          CALL DSCAL( N, TSCAL, CNORM, 1 ) >*/
        dscal_(n, &tscal, &cnorm[1], &c__1);
/*<       END IF >*/
    }

/*     Compute a bound on the computed solution vector to see if the */
/*     Level 2 BLAS routine DTRSV can be used. */

/*<       J = IDAMAX( N, X, 1 ) >*/
    j = idamax_(n, &x[1], &c__1);
/*<       XMAX = ABS( X( J ) ) >*/
    xmax = (d__1 = x[j], abs(d__1));
/*<       XBND = XMAX >*/
    xbnd = xmax;
/*<       IF( NOTRAN ) THEN >*/
    if (notran) {

/*        Compute the growth in A * x = b. */

/*<          IF( UPPER ) THEN >*/
        if (upper) {
/*<             JFIRST = N >*/
            jfirst = *n;
/*<             JLAST = 1 >*/
            jlast = 1;
/*<             JINC = -1 >*/
            jinc = -1;
/*<          ELSE >*/
        } else {
/*<             JFIRST = 1 >*/
            jfirst = 1;
/*<             JLAST = N >*/
            jlast = *n;
/*<             JINC = 1 >*/
            jinc = 1;
/*<          END IF >*/
        }

/*<          IF( TSCAL.NE.ONE ) THEN >*/
        if (tscal != 1.) {
/*<             GROW = ZERO >*/
            grow = 0.;
/*<             GO TO 50 >*/
            goto L50;
/*<          END IF >*/
        }

/*<          IF( NOUNIT ) THEN >*/
        if (nounit) {

/*           A is non-unit triangular. */

/*           Compute GROW = 1/G(j) and XBND = 1/M(j). */
/*           Initially, G(0) = max{x(i), i=1,...,n}. */

/*<             GROW = ONE / MAX( XBND, SMLNUM ) >*/
            grow = 1. / max(xbnd,smlnum);
/*<             XBND = GROW >*/
            xbnd = grow;
/*<             DO 30 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) {

/*              Exit the loop if the growth factor is too small. */

/*<    >*/
                if (grow <= smlnum) {
                    goto L50;
                }

/*              M(j) = G(j-1) / abs(A(j,j)) */

/*<                TJJ = ABS( A( J, J ) ) >*/
                tjj = (d__1 = a[j + j * a_dim1], abs(d__1));
/*<                XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) >*/
/* Computing MIN */
                d__1 = xbnd, d__2 = min(1.,tjj) * grow;
                xbnd = min(d__1,d__2);
/*<                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN >*/
                if (tjj + cnorm[j] >= smlnum) {

/*                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) */

/*<                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) >*/
                    grow *= tjj / (tjj + cnorm[j]);
/*<                ELSE >*/
                } else {

/*                 G(j) could overflow, set GROW to 0. */

/*<                   GROW = ZERO >*/
                    grow = 0.;
/*<                END IF >*/
                }
/*<    30       CONTINUE >*/
/* L30: */
            }
/*<             GROW = XBND >*/
            grow = xbnd;
/*<          ELSE >*/
        } else {

/*           A is unit triangular. */

/*           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */

/*<             GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) >*/
/* Computing MIN */
            d__1 = 1., d__2 = 1. / max(xbnd,smlnum);
            grow = min(d__1,d__2);
/*<             DO 40 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) {

/*              Exit the loop if the growth factor is too small. */

/*<    >*/
                if (grow <= smlnum) {
                    goto L50;
                }

/*              G(j) = G(j-1)*( 1 + CNORM(j) ) */

/*<                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) >*/
                grow *= 1. / (cnorm[j] + 1.);
/*<    40       CONTINUE >*/
/* L40: */
            }
/*<          END IF >*/
        }
/*<    50    CONTINUE >*/
L50:

/*<       ELSE >*/
        ;
    } else {

/*        Compute the growth in A' * x = b. */

/*<          IF( UPPER ) THEN >*/
        if (upper) {
/*<             JFIRST = 1 >*/
            jfirst = 1;
/*<             JLAST = N >*/
            jlast = *n;
/*<             JINC = 1 >*/
            jinc = 1;
/*<          ELSE >*/
        } else {
/*<             JFIRST = N >*/
            jfirst = *n;
/*<             JLAST = 1 >*/
            jlast = 1;
/*<             JINC = -1 >*/
            jinc = -1;
/*<          END IF >*/
        }

/*<          IF( TSCAL.NE.ONE ) THEN >*/
        if (tscal != 1.) {
/*<             GROW = ZERO >*/
            grow = 0.;
/*<             GO TO 80 >*/
            goto L80;
/*<          END IF >*/
        }

/*<          IF( NOUNIT ) THEN >*/
        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 = ONE / MAX( XBND, SMLNUM ) >*/
            grow = 1. / max(xbnd,smlnum);
/*<             XBND = GROW >*/
            xbnd = grow;
/*<             DO 60 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) {

/*              Exit the loop if the growth factor is too small. */

/*<    >*/
                if (grow <= smlnum) {
                    goto L80;
                }

/*              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) */

/*<                XJ = ONE + CNORM( J ) >*/
                xj = cnorm[j] + 1.;
/*<                GROW = MIN( GROW, XBND / XJ ) >*/
/* Computing MIN */
                d__1 = grow, d__2 = xbnd / xj;
                grow = min(d__1,d__2);

/*              M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) */

/*<                TJJ = ABS( A( J, J ) ) >*/
                tjj = (d__1 = a[j + j * a_dim1], abs(d__1));
/*<    >*/
                if (xj > tjj) {
                    xbnd *= tjj / xj;
                }
/*<    60       CONTINUE >*/
/* L60: */
            }
/*<             GROW = MIN( GROW, XBND ) >*/
            grow = min(grow,xbnd);
/*<          ELSE >*/
        } else {

/*           A is unit triangular. */

/*           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */

/*<             GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) >*/
/* Computing MIN */
            d__1 = 1., d__2 = 1. / max(xbnd,smlnum);
            grow = min(d__1,d__2);
/*<             DO 70 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) {

/*              Exit the loop if the growth factor is too small. */

/*<    >*/
                if (grow <= smlnum) {
                    goto L80;
                }

/*              G(j) = ( 1 + CNORM(j) )*G(j-1) */

/*<                XJ = ONE + CNORM( J ) >*/
                xj = cnorm[j] + 1.;
/*<                GROW = GROW / XJ >*/
                grow /= xj;
/*<    70       CONTINUE >*/
/* L70: */
            }
/*<          END IF >*/
        }
/*<    80    CONTINUE >*/
L80:
/*<       END IF >*/
        ;
    }

/*<       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN >*/
    if (grow * tscal > smlnum) {

/*        Use the Level 2 BLAS solve if the reciprocal of the bound on */
/*        elements of X is not too small. */

/*<          CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 ) >*/
        dtrsv_(uplo, trans, diag, n, &a[a_offset], lda, &x[1], &c__1, (ftnlen)
                1, (ftnlen)1, (ftnlen)1);
/*<       ELSE >*/
    } else {

/*        Use a Level 1 BLAS solve, scaling intermediate results. */

/*<          IF( XMAX.GT.BIGNUM ) THEN >*/
        if (xmax > bignum) {

/*           Scale X so that its components are less than or equal to */
/*           BIGNUM in absolute value. */

/*<             SCALE = BIGNUM / XMAX >*/
            *scale = bignum / xmax;
/*<             CALL DSCAL( N, SCALE, X, 1 ) >*/
            dscal_(n, scale, &x[1], &c__1);
/*<             XMAX = BIGNUM >*/
            xmax = bignum;
/*<          END IF >*/
        }

/*<          IF( NOTRAN ) THEN >*/
        if (notran) {

/*           Solve A * x = b */

/*<             DO 110 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) / A(j,j), scaling x if necessary. */

/*<                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 L100;
                    }
/*<                END IF >*/
                }
/*<                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/b(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;
/*<                   XJ = ABS( X( J ) ) >*/
                    xj = (d__1 = x[j], abs(d__1));
/*<                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 DSCAL( N, REC, X, 1 ) >*/
                        dscal_(n, &rec, &x[1], &c__1);
/*<                      SCALE = SCALE*REC >*/
                        *scale *= rec;
/*<                      XMAX = XMAX*REC >*/