zlatrs.c

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

            }
/*<             CNORM( N ) = ZERO >*/
            cnorm[*n] = 0.;
/*<          END IF >*/
        }
/*<       END IF >*/
    }

/*     Scale the column norms by TSCAL if the maximum element in CNORM is */
/*     greater than BIGNUM/2. */

/*<       IMAX = IDAMAX( N, CNORM, 1 ) >*/
    imax = idamax_(n, &cnorm[1], &c__1);
/*<       TMAX = CNORM( IMAX ) >*/
    tmax = cnorm[imax];
/*<       IF( TMAX.LE.BIGNUM*HALF ) THEN >*/
    if (tmax <= bignum * .5) {
/*<          TSCAL = ONE >*/
        tscal = 1.;
/*<       ELSE >*/
    } else {
/*<          TSCAL = HALF / ( SMLNUM*TMAX ) >*/
        tscal = .5 / (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 ZTRSV can be used. */

/*<       XMAX = ZERO >*/
    xmax = 0.;
/*<       DO 30 J = 1, N >*/
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/*<          XMAX = MAX( XMAX, CABS2( X( J ) ) ) >*/
/* Computing MAX */
        i__2 = j;
        d__3 = xmax, d__4 = (d__1 = x[i__2].r / 2., abs(d__1)) + (d__2 = 
                d_imag(&x[j]) / 2., abs(d__2));
        xmax = max(d__3,d__4);
/*<    30 CONTINUE >*/
/* L30: */
    }
/*<       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 60 >*/
            goto L60;
/*<          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 = HALF / MAX( XBND, SMLNUM ) >*/
            grow = .5 / max(xbnd,smlnum);
/*<             XBND = GROW >*/
            xbnd = grow;
/*<             DO 40 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 L60;
                }

/*<                TJJS = A( J, J ) >*/
                i__3 = j + j * a_dim1;
                tjjs.r = a[i__3].r, tjjs.i = a[i__3].i;
/*<                TJJ = CABS1( TJJS ) >*/
                tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(
                        d__2));

/*<                IF( TJJ.GE.SMLNUM ) THEN >*/
                if (tjj >= smlnum) {

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

/*<                   XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) >*/
/* Computing MIN */
                    d__1 = xbnd, d__2 = min(1.,tjj) * grow;
                    xbnd = min(d__1,d__2);
/*<                ELSE >*/
                } else {

/*                 M(j) could overflow, set XBND to 0. */

/*<                   XBND = ZERO >*/
                    xbnd = 0.;
/*<                END IF >*/
                }

/*<                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 >*/
                }
/*<    40       CONTINUE >*/
/* L40: */
            }
/*<             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, HALF / MAX( XBND, SMLNUM ) ) >*/
/* Computing MIN */
            d__1 = 1., d__2 = .5 / max(xbnd,smlnum);
            grow = min(d__1,d__2);
/*<             DO 50 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 L60;
                }

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

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

/*<       ELSE >*/
        ;
    } else {

/*        Compute the growth in A**T * x = b  or  A**H * 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 90 >*/
            goto L90;
/*<          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 = HALF / MAX( XBND, SMLNUM ) >*/
            grow = .5 / max(xbnd,smlnum);
/*<             XBND = GROW >*/
            xbnd = grow;
/*<             DO 70 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 L90;
                }

/*              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);

/*<                TJJS = A( J, J ) >*/
                i__3 = j + j * a_dim1;
                tjjs.r = a[i__3].r, tjjs.i = a[i__3].i;
/*<                TJJ = CABS1( TJJS ) >*/
                tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(
                        d__2));

/*<                IF( TJJ.GE.SMLNUM ) THEN >*/
                if (tjj >= smlnum) {

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

/*<    >*/
                    if (xj > tjj) {
                        xbnd *= tjj / xj;
                    }
/*<                ELSE >*/
                } else {

/*                 M(j) could overflow, set XBND to 0. */

/*<                   XBND = ZERO >*/
                    xbnd = 0.;
/*<                END IF >*/
                }
/*<    70       CONTINUE >*/
/* L70: */
            }
/*<             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, HALF / MAX( XBND, SMLNUM ) ) >*/
/* Computing MIN */
            d__1 = 1., d__2 = .5 / max(xbnd,smlnum);
            grow = min(d__1,d__2);
/*<             DO 80 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 L90;
                }

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

/*<                XJ = ONE + CNORM( J ) >*/
                xj = cnorm[j] + 1.;
/*<                GROW = GROW / XJ >*/
                grow /= xj;
/*<    80       CONTINUE >*/
/* L80: */
            }
/*<          END IF >*/
        }
/*<    90    CONTINUE >*/
L90:
/*<       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 ZTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 ) >*/
        ztrsv_(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*HALF ) THEN >*/
        if (xmax > bignum * .5) {

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

/*<             SCALE = ( BIGNUM*HALF ) / XMAX >*/
            *scale = bignum * .5 / xmax;
/*<             CALL ZDSCAL( N, SCALE, X, 1 ) >*/
            zdscal_(n, scale, &x[1], &c__1);
/*<             XMAX = BIGNUM >*/
            xmax = bignum;
/*<          ELSE >*/
        } else {
/*<             XMAX = XMAX*TWO >*/
            xmax *= 2.;
/*<          END IF >*/
        }

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

/*           Solve A * x = b */

/*<             DO 120 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 = 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 = A( J, J )*TSCAL >*/
                    i__3 = j + j * a_dim1;
                    z__1.r = tscal * a[i__3].r, z__1.i = tscal * a[i__3].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 L110;
                    }
/*<                END IF >*/
                }
/*<                TJJ = CABS1( TJJS ) >*/
                tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(
                        d__2));
/*<                IF( TJJ.GT.SMLNUM ) THEN >*/