diff options
author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
---|---|---|
committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/zlatbs.f | |
download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip |
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/zlatbs.f')
-rw-r--r-- | SRC/zlatbs.f | 908 |
1 files changed, 908 insertions, 0 deletions
diff --git a/SRC/zlatbs.f b/SRC/zlatbs.f new file mode 100644 index 00000000..5dccd835 --- /dev/null +++ b/SRC/zlatbs.f @@ -0,0 +1,908 @@ + SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, + $ SCALE, CNORM, INFO ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER DIAG, NORMIN, TRANS, UPLO + INTEGER INFO, KD, LDAB, N + DOUBLE PRECISION SCALE +* .. +* .. Array Arguments .. + DOUBLE PRECISION CNORM( * ) + COMPLEX*16 AB( LDAB, * ), X( * ) +* .. +* +* Purpose +* ======= +* +* ZLATBS solves one of the triangular systems +* +* A * x = s*b, A**T * x = s*b, or A**H * x = s*b, +* +* with scaling to prevent overflow, where A is an upper or lower +* triangular band matrix. Here A' denotes the transpose of A, x and b +* are n-element vectors, and s is a scaling factor, usually less than +* or equal to 1, chosen so that the components of x will be less than +* the overflow threshold. If the unscaled problem will not cause +* overflow, the Level 2 BLAS routine ZTBSV is called. If the matrix A +* is singular (A(j,j) = 0 for some j), then s is set to 0 and a +* non-trivial solution to A*x = 0 is returned. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the matrix A is upper or lower triangular. +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* TRANS (input) CHARACTER*1 +* Specifies the operation applied to A. +* = 'N': Solve A * x = s*b (No transpose) +* = 'T': Solve A**T * x = s*b (Transpose) +* = 'C': Solve A**H * x = s*b (Conjugate transpose) +* +* DIAG (input) CHARACTER*1 +* Specifies whether or not the matrix A is unit triangular. +* = 'N': Non-unit triangular +* = 'U': Unit triangular +* +* NORMIN (input) CHARACTER*1 +* Specifies whether CNORM has been set or not. +* = 'Y': CNORM contains the column norms on entry +* = 'N': CNORM is not set on entry. On exit, the norms will +* be computed and stored in CNORM. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* KD (input) INTEGER +* The number of subdiagonals or superdiagonals in the +* triangular matrix A. KD >= 0. +* +* AB (input) COMPLEX*16 array, dimension (LDAB,N) +* The upper or lower triangular band matrix A, stored in the +* first KD+1 rows of the array. The j-th column of A is stored +* in the j-th column of the array AB as follows: +* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; +* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). +* +* LDAB (input) INTEGER +* The leading dimension of the array AB. LDAB >= KD+1. +* +* X (input/output) COMPLEX*16 array, dimension (N) +* On entry, the right hand side b of the triangular system. +* On exit, X is overwritten by the solution vector x. +* +* SCALE (output) DOUBLE PRECISION +* The scaling factor s for the triangular system +* A * x = s*b, A**T * x = s*b, or A**H * x = s*b. +* If SCALE = 0, the matrix A is singular or badly scaled, and +* the vector x is an exact or approximate solution to A*x = 0. +* +* CNORM (input or output) DOUBLE PRECISION array, dimension (N) +* +* If NORMIN = 'Y', CNORM is an input argument and CNORM(j) +* contains the norm of the off-diagonal part of the j-th column +* of A. If TRANS = 'N', CNORM(j) must be greater than or equal +* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) +* must be greater than or equal to the 1-norm. +* +* If NORMIN = 'N', CNORM is an output argument and CNORM(j) +* returns the 1-norm of the offdiagonal part of the j-th column +* of A. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -k, the k-th argument had an illegal value +* +* Further Details +* ======= ======= +* +* A rough bound on x is computed; if that is less than overflow, ZTBSV +* is called, otherwise, specific code is used which checks for possible +* overflow or divide-by-zero at every operation. +* +* A columnwise scheme is used for solving A*x = b. The basic algorithm +* if A is lower triangular is +* +* x[1:n] := b[1:n] +* for j = 1, ..., n +* x(j) := x(j) / A(j,j) +* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] +* end +* +* Define bounds on the components of x after j iterations of the loop: +* M(j) = bound on x[1:j] +* G(j) = bound on x[j+1:n] +* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. +* +* Then for iteration j+1 we have +* M(j+1) <= G(j) / | A(j+1,j+1) | +* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | +* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) +* +* where CNORM(j+1) is greater than or equal to the infinity-norm of +* column j+1 of A, not counting the diagonal. Hence +* +* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) +* 1<=i<=j +* and +* +* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) +* 1<=i< j +* +* Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTBSV if the +* reciprocal of the largest M(j), j=1,..,n, is larger than +* max(underflow, 1/overflow). +* +* The bound on x(j) is also used to determine when a step in the +* columnwise method can be performed without fear of overflow. If +* the computed bound is greater than a large constant, x is scaled to +* prevent overflow, but if the bound overflows, x is set to 0, x(j) to +* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. +* +* Similarly, a row-wise scheme is used to solve A**T *x = b or +* A**H *x = b. The basic algorithm for A upper triangular is +* +* for j = 1, ..., n +* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) +* end +* +* We simultaneously compute two bounds +* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j +* M(j) = bound on x(i), 1<=i<=j +* +* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we +* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. +* Then the bound on x(j) is +* +* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | +* +* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) +* 1<=i<=j +* +* and we can safely call ZTBSV if 1/M(n) and 1/G(n) are both greater +* than max(underflow, 1/overflow). +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, HALF, ONE, TWO + PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0, + $ TWO = 2.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL NOTRAN, NOUNIT, UPPER + INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND + DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL, + $ XBND, XJ, XMAX + COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER IDAMAX, IZAMAX + DOUBLE PRECISION DLAMCH, DZASUM + COMPLEX*16 ZDOTC, ZDOTU, ZLADIV + EXTERNAL LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC, + $ ZDOTU, ZLADIV +* .. +* .. External Subroutines .. + EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN +* .. +* .. Statement Functions .. + DOUBLE PRECISION CABS1, CABS2 +* .. +* .. Statement Function definitions .. + CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) + CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) + + $ ABS( DIMAG( ZDUM ) / 2.D0 ) +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + NOTRAN = LSAME( TRANS, 'N' ) + NOUNIT = LSAME( DIAG, 'N' ) +* +* Test the input parameters. +* + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. + $ LSAME( TRANS, 'C' ) ) THEN + INFO = -2 + ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN + INFO = -3 + ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. + $ LSAME( NORMIN, 'N' ) ) THEN + INFO = -4 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( KD.LT.0 ) THEN + INFO = -6 + ELSE IF( LDAB.LT.KD+1 ) THEN + INFO = -8 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZLATBS', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine machine dependent parameters to control overflow. +* + SMLNUM = DLAMCH( 'Safe minimum' ) + BIGNUM = ONE / SMLNUM + CALL DLABAD( SMLNUM, BIGNUM ) + SMLNUM = SMLNUM / DLAMCH( 'Precision' ) + BIGNUM = ONE / SMLNUM + SCALE = ONE +* + IF( LSAME( NORMIN, 'N' ) ) THEN +* +* Compute the 1-norm of each column, not including the diagonal. +* + IF( UPPER ) THEN +* +* A is upper triangular. +* + DO 10 J = 1, N + JLEN = MIN( KD, J-1 ) + CNORM( J ) = DZASUM( JLEN, AB( KD+1-JLEN, J ), 1 ) + 10 CONTINUE + ELSE +* +* A is lower triangular. +* + DO 20 J = 1, N + JLEN = MIN( KD, N-J ) + IF( JLEN.GT.0 ) THEN + CNORM( J ) = DZASUM( JLEN, AB( 2, J ), 1 ) + ELSE + CNORM( J ) = ZERO + END IF + 20 CONTINUE + 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 ) + TMAX = CNORM( IMAX ) + IF( TMAX.LE.BIGNUM*HALF ) THEN + TSCAL = ONE + ELSE + TSCAL = HALF / ( SMLNUM*TMAX ) + CALL DSCAL( N, TSCAL, CNORM, 1 ) + END IF +* +* Compute a bound on the computed solution vector to see if the +* Level 2 BLAS routine ZTBSV can be used. +* + XMAX = ZERO + DO 30 J = 1, N + XMAX = MAX( XMAX, CABS2( X( J ) ) ) + 30 CONTINUE + XBND = XMAX + IF( NOTRAN ) THEN +* +* Compute the growth in A * x = b. +* + IF( UPPER ) THEN + JFIRST = N + JLAST = 1 + JINC = -1 + MAIND = KD + 1 + ELSE + JFIRST = 1 + JLAST = N + JINC = 1 + MAIND = 1 + END IF +* + IF( TSCAL.NE.ONE ) THEN + GROW = ZERO + GO TO 60 + END IF +* + IF( NOUNIT ) THEN +* +* 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 ) + XBND = GROW + DO 40 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 60 +* + TJJS = AB( MAIND, J ) + TJJ = CABS1( TJJS ) +* + IF( TJJ.GE.SMLNUM ) THEN +* +* M(j) = G(j-1) / abs(A(j,j)) +* + XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) + ELSE +* +* M(j) could overflow, set XBND to 0. +* + XBND = ZERO + END IF +* + IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN +* +* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) +* + GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) + ELSE +* +* G(j) could overflow, set GROW to 0. +* + GROW = ZERO + END IF + 40 CONTINUE + GROW = XBND + 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 ) ) + DO 50 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 60 +* +* G(j) = G(j-1)*( 1 + CNORM(j) ) +* + GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) + 50 CONTINUE + END IF + 60 CONTINUE +* + ELSE +* +* Compute the growth in A**T * x = b or A**H * x = b. +* + IF( UPPER ) THEN + JFIRST = 1 + JLAST = N + JINC = 1 + MAIND = KD + 1 + ELSE + JFIRST = N + JLAST = 1 + JINC = -1 + MAIND = 1 + END IF +* + IF( TSCAL.NE.ONE ) THEN + GROW = ZERO + GO TO 90 + END IF +* + IF( NOUNIT ) THEN +* +* 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 ) + XBND = GROW + DO 70 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 90 +* +* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) +* + XJ = ONE + CNORM( J ) + GROW = MIN( GROW, XBND / XJ ) +* + TJJS = AB( MAIND, J ) + TJJ = CABS1( TJJS ) +* + IF( TJJ.GE.SMLNUM ) THEN +* +* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) +* + IF( XJ.GT.TJJ ) + $ XBND = XBND*( TJJ / XJ ) + ELSE +* +* M(j) could overflow, set XBND to 0. +* + XBND = ZERO + END IF + 70 CONTINUE + GROW = MIN( GROW, XBND ) + 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 ) ) + DO 80 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 90 +* +* G(j) = ( 1 + CNORM(j) )*G(j-1) +* + XJ = ONE + CNORM( J ) + GROW = GROW / XJ + 80 CONTINUE + END IF + 90 CONTINUE + END IF +* + IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN +* +* Use the Level 2 BLAS solve if the reciprocal of the bound on +* elements of X is not too small. +* + CALL ZTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 ) + ELSE +* +* Use a Level 1 BLAS solve, scaling intermediate results. +* + IF( XMAX.GT.BIGNUM*HALF ) THEN +* +* Scale X so that its components are less than or equal to +* BIGNUM in absolute value. +* + SCALE = ( BIGNUM*HALF ) / XMAX + CALL ZDSCAL( N, SCALE, X, 1 ) + XMAX = BIGNUM + ELSE + XMAX = XMAX*TWO + END IF +* + IF( NOTRAN ) THEN +* +* Solve A * x = b +* + DO 120 J = JFIRST, JLAST, JINC +* +* Compute x(j) = b(j) / A(j,j), scaling x if necessary. +* + XJ = CABS1( X( J ) ) + IF( NOUNIT ) THEN + TJJS = AB( MAIND, J )*TSCAL + ELSE + TJJS = TSCAL + IF( TSCAL.EQ.ONE ) + $ GO TO 110 + END IF + TJJ = CABS1( TJJS ) + IF( TJJ.GT.SMLNUM ) THEN +* +* abs(A(j,j)) > SMLNUM: +* + IF( TJJ.LT.ONE ) THEN + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by 1/b(j). +* + REC = ONE / XJ + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF + X( J ) = ZLADIV( X( J ), TJJS ) + XJ = CABS1( X( J ) ) + ELSE IF( TJJ.GT.ZERO ) THEN +* +* 0 < abs(A(j,j)) <= SMLNUM: +* + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* 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 + IF( CNORM( J ).GT.ONE ) THEN +* +* Scale by 1/CNORM(j) to avoid overflow when +* multiplying x(j) times column j. +* + REC = REC / CNORM( J ) + END IF + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + X( J ) = ZLADIV( X( J ), TJJS ) + XJ = CABS1( X( J ) ) + 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 100 I = 1, N + X( I ) = ZERO + 100 CONTINUE + X( J ) = ONE + XJ = ONE + SCALE = ZERO + XMAX = ZERO + END IF + 110 CONTINUE +* +* Scale x if necessary to avoid overflow when adding a +* multiple of column j of A. +* + IF( XJ.GT.ONE ) THEN + REC = ONE / XJ + IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN +* +* Scale x by 1/(2*abs(x(j))). +* + REC = REC*HALF + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + END IF + ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN +* +* Scale x by 1/2. +* + CALL ZDSCAL( N, HALF, X, 1 ) + SCALE = SCALE*HALF + END IF +* + IF( UPPER ) THEN + IF( J.GT.1 ) THEN +* +* Compute the update +* x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) - +* x(j)* A(max(1,j-kd):j-1,j) +* + JLEN = MIN( KD, J-1 ) + CALL ZAXPY( JLEN, -X( J )*TSCAL, + $ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 ) + I = IZAMAX( J-1, X, 1 ) + XMAX = CABS1( X( I ) ) + END IF + ELSE IF( J.LT.N ) THEN +* +* Compute the update +* x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) - +* x(j) * A(j+1:min(j+kd,n),j) +* + JLEN = MIN( KD, N-J ) + IF( JLEN.GT.0 ) + $ CALL ZAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1, + $ X( J+1 ), 1 ) + I = J + IZAMAX( N-J, X( J+1 ), 1 ) + XMAX = CABS1( X( I ) ) + END IF + 120 CONTINUE +* + ELSE IF( LSAME( TRANS, 'T' ) ) THEN +* +* Solve A**T * x = b +* + DO 170 J = JFIRST, JLAST, JINC +* +* Compute x(j) = b(j) - sum A(k,j)*x(k). +* k<>j +* + XJ = CABS1( X( J ) ) + USCAL = TSCAL + REC = ONE / MAX( XMAX, ONE ) + IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN +* +* If x(j) could overflow, scale x by 1/(2*XMAX). +* + REC = REC*HALF + IF( NOUNIT ) THEN + TJJS = AB( MAIND, J )*TSCAL + ELSE + TJJS = TSCAL + END IF + TJJ = CABS1( TJJS ) + IF( TJJ.GT.ONE ) THEN +* +* Divide by A(j,j) when scaling x if A(j,j) > 1. +* + REC = MIN( ONE, REC*TJJ ) + USCAL = ZLADIV( USCAL, TJJS ) + END IF + IF( REC.LT.ONE ) THEN + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF +* + CSUMJ = ZERO + IF( USCAL.EQ.DCMPLX( ONE ) ) THEN +* +* If the scaling needed for A in the dot product is 1, +* call ZDOTU to perform the dot product. +* + IF( UPPER ) THEN + JLEN = MIN( KD, J-1 ) + CSUMJ = ZDOTU( JLEN, AB( KD+1-JLEN, J ), 1, + $ X( J-JLEN ), 1 ) + ELSE + JLEN = MIN( KD, N-J ) + IF( JLEN.GT.1 ) + $ CSUMJ = ZDOTU( JLEN, AB( 2, J ), 1, X( J+1 ), + $ 1 ) + END IF + ELSE +* +* Otherwise, use in-line code for the dot product. +* + IF( UPPER ) THEN + JLEN = MIN( KD, J-1 ) + DO 130 I = 1, JLEN + CSUMJ = CSUMJ + ( AB( KD+I-JLEN, J )*USCAL )* + $ X( J-JLEN-1+I ) + 130 CONTINUE + ELSE + JLEN = MIN( KD, N-J ) + DO 140 I = 1, JLEN + CSUMJ = CSUMJ + ( AB( I+1, J )*USCAL )*X( J+I ) + 140 CONTINUE + END IF + END IF +* + IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN +* +* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) +* was not used to scale the dotproduct. +* + X( J ) = X( J ) - CSUMJ + XJ = CABS1( X( J ) ) + IF( NOUNIT ) THEN +* +* Compute x(j) = x(j) / A(j,j), scaling if necessary. +* + TJJS = AB( MAIND, J )*TSCAL + ELSE + TJJS = TSCAL + IF( TSCAL.EQ.ONE ) + $ GO TO 160 + END IF + TJJ = CABS1( TJJS ) + IF( TJJ.GT.SMLNUM ) THEN +* +* abs(A(j,j)) > SMLNUM: +* + IF( TJJ.LT.ONE ) THEN + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale X by 1/abs(x(j)). +* + REC = ONE / XJ + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF + X( J ) = ZLADIV( X( J ), TJJS ) + ELSE IF( TJJ.GT.ZERO ) THEN +* +* 0 < abs(A(j,j)) <= SMLNUM: +* + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. +* + REC = ( TJJ*BIGNUM ) / XJ + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + X( J ) = ZLADIV( 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. +* + DO 150 I = 1, N + X( I ) = ZERO + 150 CONTINUE + X( J ) = ONE + SCALE = ZERO + XMAX = ZERO + END IF + 160 CONTINUE + ELSE +* +* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot +* product has already been divided by 1/A(j,j). +* + X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ + END IF + XMAX = MAX( XMAX, CABS1( X( J ) ) ) + 170 CONTINUE +* + ELSE +* +* Solve A**H * x = b +* + DO 220 J = JFIRST, JLAST, JINC +* +* Compute x(j) = b(j) - sum A(k,j)*x(k). +* k<>j +* + XJ = CABS1( X( J ) ) + USCAL = TSCAL + REC = ONE / MAX( XMAX, ONE ) + IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN +* +* If x(j) could overflow, scale x by 1/(2*XMAX). +* + REC = REC*HALF + IF( NOUNIT ) THEN + TJJS = DCONJG( AB( MAIND, J ) )*TSCAL + ELSE + TJJS = TSCAL + END IF + TJJ = CABS1( TJJS ) + IF( TJJ.GT.ONE ) THEN +* +* Divide by A(j,j) when scaling x if A(j,j) > 1. +* + REC = MIN( ONE, REC*TJJ ) + USCAL = ZLADIV( USCAL, TJJS ) + END IF + IF( REC.LT.ONE ) THEN + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF +* + CSUMJ = ZERO + IF( USCAL.EQ.DCMPLX( ONE ) ) THEN +* +* If the scaling needed for A in the dot product is 1, +* call ZDOTC to perform the dot product. +* + IF( UPPER ) THEN + JLEN = MIN( KD, J-1 ) + CSUMJ = ZDOTC( JLEN, AB( KD+1-JLEN, J ), 1, + $ X( J-JLEN ), 1 ) + ELSE + JLEN = MIN( KD, N-J ) + IF( JLEN.GT.1 ) + $ CSUMJ = ZDOTC( JLEN, AB( 2, J ), 1, X( J+1 ), + $ 1 ) + END IF + ELSE +* +* Otherwise, use in-line code for the dot product. +* + IF( UPPER ) THEN + JLEN = MIN( KD, J-1 ) + DO 180 I = 1, JLEN + CSUMJ = CSUMJ + ( DCONJG( AB( KD+I-JLEN, J ) )* + $ USCAL )*X( J-JLEN-1+I ) + 180 CONTINUE + ELSE + JLEN = MIN( KD, N-J ) + DO 190 I = 1, JLEN + CSUMJ = CSUMJ + ( DCONJG( AB( I+1, J ) )*USCAL ) + $ *X( J+I ) + 190 CONTINUE + END IF + END IF +* + IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN +* +* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) +* was not used to scale the dotproduct. +* + X( J ) = X( J ) - CSUMJ + XJ = CABS1( X( J ) ) + IF( NOUNIT ) THEN +* +* Compute x(j) = x(j) / A(j,j), scaling if necessary. +* + TJJS = DCONJG( AB( MAIND, J ) )*TSCAL + ELSE + TJJS = TSCAL + IF( TSCAL.EQ.ONE ) + $ GO TO 210 + END IF + TJJ = CABS1( TJJS ) + IF( TJJ.GT.SMLNUM ) THEN +* +* abs(A(j,j)) > SMLNUM: +* + IF( TJJ.LT.ONE ) THEN + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale X by 1/abs(x(j)). +* + REC = ONE / XJ + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF + X( J ) = ZLADIV( X( J ), TJJS ) + ELSE IF( TJJ.GT.ZERO ) THEN +* +* 0 < abs(A(j,j)) <= SMLNUM: +* + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. +* + REC = ( TJJ*BIGNUM ) / XJ + CALL ZDSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + X( J ) = ZLADIV( 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. +* + DO 200 I = 1, N + X( I ) = ZERO + 200 CONTINUE + X( J ) = ONE + SCALE = ZERO + XMAX = ZERO + END IF + 210 CONTINUE + ELSE +* +* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot +* product has already been divided by 1/A(j,j). +* + X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ + END IF + XMAX = MAX( XMAX, CABS1( X( J ) ) ) + 220 CONTINUE + END IF + SCALE = SCALE / TSCAL + END IF +* +* Scale the column norms by 1/TSCAL for return. +* + IF( TSCAL.NE.ONE ) THEN + CALL DSCAL( N, ONE / TSCAL, CNORM, 1 ) + END IF +* + RETURN +* +* End of ZLATBS +* + END |