LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zlatbs()

 subroutine zlatbs ( character uplo, character trans, character diag, character normin, integer n, integer kd, complex*16, dimension( ldab, * ) ab, integer ldab, complex*16, dimension( * ) x, double precision scale, double precision, dimension( * ) cnorm, integer info )

ZLATBS solves a triangular banded system of equations.

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**T 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.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangular``` [in] TRANS ``` TRANS is 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)``` [in] DIAG ``` DIAG is CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangular``` [in] NORMIN ``` NORMIN is 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.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] KD ``` KD is INTEGER The number of subdiagonals or superdiagonals in the triangular matrix A. KD >= 0.``` [in] AB ``` AB is 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).``` [in] LDAB ``` LDAB is INTEGER The leading dimension of the array AB. LDAB >= KD+1.``` [in,out] X ``` X is 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.``` [out] SCALE ``` SCALE is 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.``` [in,out] CNORM ``` CNORM is 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.``` [out] INFO ``` INFO is 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).```

Definition at line 241 of file zlatbs.f.

243*
244* -- LAPACK auxiliary routine --
245* -- LAPACK is a software package provided by Univ. of Tennessee, --
246* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
247*
248* .. Scalar Arguments ..
249 CHARACTER DIAG, NORMIN, TRANS, UPLO
250 INTEGER INFO, KD, LDAB, N
251 DOUBLE PRECISION SCALE
252* ..
253* .. Array Arguments ..
254 DOUBLE PRECISION CNORM( * )
255 COMPLEX*16 AB( LDAB, * ), X( * )
256* ..
257*
258* =====================================================================
259*
260* .. Parameters ..
261 DOUBLE PRECISION ZERO, HALF, ONE, TWO
262 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0,
263 \$ two = 2.0d+0 )
264* ..
265* .. Local Scalars ..
266 LOGICAL NOTRAN, NOUNIT, UPPER
267 INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
268 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
269 \$ XBND, XJ, XMAX
270 COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
271* ..
272* .. External Functions ..
273 LOGICAL LSAME
274 INTEGER IDAMAX, IZAMAX
275 DOUBLE PRECISION DLAMCH, DZASUM
277 EXTERNAL lsame, idamax, izamax, dlamch, dzasum, zdotc,
279* ..
280* .. External Subroutines ..
281 EXTERNAL dscal, xerbla, zaxpy, zdscal, ztbsv
282* ..
283* .. Intrinsic Functions ..
284 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
285* ..
286* .. Statement Functions ..
287 DOUBLE PRECISION CABS1, CABS2
288* ..
289* .. Statement Function definitions ..
290 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
291 cabs2( zdum ) = abs( dble( zdum ) / 2.d0 ) +
292 \$ abs( dimag( zdum ) / 2.d0 )
293* ..
294* .. Executable Statements ..
295*
296 info = 0
297 upper = lsame( uplo, 'U' )
298 notran = lsame( trans, 'N' )
299 nounit = lsame( diag, 'N' )
300*
301* Test the input parameters.
302*
303 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
304 info = -1
305 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
306 \$ lsame( trans, 'C' ) ) THEN
307 info = -2
308 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
309 info = -3
310 ELSE IF( .NOT.lsame( normin, 'Y' ) .AND. .NOT.
311 \$ lsame( normin, 'N' ) ) THEN
312 info = -4
313 ELSE IF( n.LT.0 ) THEN
314 info = -5
315 ELSE IF( kd.LT.0 ) THEN
316 info = -6
317 ELSE IF( ldab.LT.kd+1 ) THEN
318 info = -8
319 END IF
320 IF( info.NE.0 ) THEN
321 CALL xerbla( 'ZLATBS', -info )
322 RETURN
323 END IF
324*
325* Quick return if possible
326*
327 scale = one
328 IF( n.EQ.0 )
329 \$ RETURN
330*
331* Determine machine dependent parameters to control overflow.
332*
333 smlnum = dlamch( 'Safe minimum' ) / dlamch( 'Precision' )
334 bignum = one / smlnum
335*
336 IF( lsame( normin, 'N' ) ) THEN
337*
338* Compute the 1-norm of each column, not including the diagonal.
339*
340 IF( upper ) THEN
341*
342* A is upper triangular.
343*
344 DO 10 j = 1, n
345 jlen = min( kd, j-1 )
346 cnorm( j ) = dzasum( jlen, ab( kd+1-jlen, j ), 1 )
347 10 CONTINUE
348 ELSE
349*
350* A is lower triangular.
351*
352 DO 20 j = 1, n
353 jlen = min( kd, n-j )
354 IF( jlen.GT.0 ) THEN
355 cnorm( j ) = dzasum( jlen, ab( 2, j ), 1 )
356 ELSE
357 cnorm( j ) = zero
358 END IF
359 20 CONTINUE
360 END IF
361 END IF
362*
363* Scale the column norms by TSCAL if the maximum element in CNORM is
364* greater than BIGNUM/2.
365*
366 imax = idamax( n, cnorm, 1 )
367 tmax = cnorm( imax )
368 IF( tmax.LE.bignum*half ) THEN
369 tscal = one
370 ELSE
371 tscal = half / ( smlnum*tmax )
372 CALL dscal( n, tscal, cnorm, 1 )
373 END IF
374*
375* Compute a bound on the computed solution vector to see if the
376* Level 2 BLAS routine ZTBSV can be used.
377*
378 xmax = zero
379 DO 30 j = 1, n
380 xmax = max( xmax, cabs2( x( j ) ) )
381 30 CONTINUE
382 xbnd = xmax
383 IF( notran ) THEN
384*
385* Compute the growth in A * x = b.
386*
387 IF( upper ) THEN
388 jfirst = n
389 jlast = 1
390 jinc = -1
391 maind = kd + 1
392 ELSE
393 jfirst = 1
394 jlast = n
395 jinc = 1
396 maind = 1
397 END IF
398*
399 IF( tscal.NE.one ) THEN
400 grow = zero
401 GO TO 60
402 END IF
403*
404 IF( nounit ) THEN
405*
406* A is non-unit triangular.
407*
408* Compute GROW = 1/G(j) and XBND = 1/M(j).
409* Initially, G(0) = max{x(i), i=1,...,n}.
410*
411 grow = half / max( xbnd, smlnum )
412 xbnd = grow
413 DO 40 j = jfirst, jlast, jinc
414*
415* Exit the loop if the growth factor is too small.
416*
417 IF( grow.LE.smlnum )
418 \$ GO TO 60
419*
420 tjjs = ab( maind, j )
421 tjj = cabs1( tjjs )
422*
423 IF( tjj.GE.smlnum ) THEN
424*
425* M(j) = G(j-1) / abs(A(j,j))
426*
427 xbnd = min( xbnd, min( one, tjj )*grow )
428 ELSE
429*
430* M(j) could overflow, set XBND to 0.
431*
432 xbnd = zero
433 END IF
434*
435 IF( tjj+cnorm( j ).GE.smlnum ) THEN
436*
437* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
438*
439 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
440 ELSE
441*
442* G(j) could overflow, set GROW to 0.
443*
444 grow = zero
445 END IF
446 40 CONTINUE
447 grow = xbnd
448 ELSE
449*
450* A is unit triangular.
451*
452* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
453*
454 grow = min( one, half / max( xbnd, smlnum ) )
455 DO 50 j = jfirst, jlast, jinc
456*
457* Exit the loop if the growth factor is too small.
458*
459 IF( grow.LE.smlnum )
460 \$ GO TO 60
461*
462* G(j) = G(j-1)*( 1 + CNORM(j) )
463*
464 grow = grow*( one / ( one+cnorm( j ) ) )
465 50 CONTINUE
466 END IF
467 60 CONTINUE
468*
469 ELSE
470*
471* Compute the growth in A**T * x = b or A**H * x = b.
472*
473 IF( upper ) THEN
474 jfirst = 1
475 jlast = n
476 jinc = 1
477 maind = kd + 1
478 ELSE
479 jfirst = n
480 jlast = 1
481 jinc = -1
482 maind = 1
483 END IF
484*
485 IF( tscal.NE.one ) THEN
486 grow = zero
487 GO TO 90
488 END IF
489*
490 IF( nounit ) THEN
491*
492* A is non-unit triangular.
493*
494* Compute GROW = 1/G(j) and XBND = 1/M(j).
495* Initially, M(0) = max{x(i), i=1,...,n}.
496*
497 grow = half / max( xbnd, smlnum )
498 xbnd = grow
499 DO 70 j = jfirst, jlast, jinc
500*
501* Exit the loop if the growth factor is too small.
502*
503 IF( grow.LE.smlnum )
504 \$ GO TO 90
505*
506* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
507*
508 xj = one + cnorm( j )
509 grow = min( grow, xbnd / xj )
510*
511 tjjs = ab( maind, j )
512 tjj = cabs1( tjjs )
513*
514 IF( tjj.GE.smlnum ) THEN
515*
516* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
517*
518 IF( xj.GT.tjj )
519 \$ xbnd = xbnd*( tjj / xj )
520 ELSE
521*
522* M(j) could overflow, set XBND to 0.
523*
524 xbnd = zero
525 END IF
526 70 CONTINUE
527 grow = min( grow, xbnd )
528 ELSE
529*
530* A is unit triangular.
531*
532* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
533*
534 grow = min( one, half / max( xbnd, smlnum ) )
535 DO 80 j = jfirst, jlast, jinc
536*
537* Exit the loop if the growth factor is too small.
538*
539 IF( grow.LE.smlnum )
540 \$ GO TO 90
541*
542* G(j) = ( 1 + CNORM(j) )*G(j-1)
543*
544 xj = one + cnorm( j )
545 grow = grow / xj
546 80 CONTINUE
547 END IF
548 90 CONTINUE
549 END IF
550*
551 IF( ( grow*tscal ).GT.smlnum ) THEN
552*
553* Use the Level 2 BLAS solve if the reciprocal of the bound on
554* elements of X is not too small.
555*
556 CALL ztbsv( uplo, trans, diag, n, kd, ab, ldab, x, 1 )
557 ELSE
558*
559* Use a Level 1 BLAS solve, scaling intermediate results.
560*
561 IF( xmax.GT.bignum*half ) THEN
562*
563* Scale X so that its components are less than or equal to
564* BIGNUM in absolute value.
565*
566 scale = ( bignum*half ) / xmax
567 CALL zdscal( n, scale, x, 1 )
568 xmax = bignum
569 ELSE
570 xmax = xmax*two
571 END IF
572*
573 IF( notran ) THEN
574*
575* Solve A * x = b
576*
577 DO 120 j = jfirst, jlast, jinc
578*
579* Compute x(j) = b(j) / A(j,j), scaling x if necessary.
580*
581 xj = cabs1( x( j ) )
582 IF( nounit ) THEN
583 tjjs = ab( maind, j )*tscal
584 ELSE
585 tjjs = tscal
586 IF( tscal.EQ.one )
587 \$ GO TO 110
588 END IF
589 tjj = cabs1( tjjs )
590 IF( tjj.GT.smlnum ) THEN
591*
592* abs(A(j,j)) > SMLNUM:
593*
594 IF( tjj.LT.one ) THEN
595 IF( xj.GT.tjj*bignum ) THEN
596*
597* Scale x by 1/b(j).
598*
599 rec = one / xj
600 CALL zdscal( n, rec, x, 1 )
601 scale = scale*rec
602 xmax = xmax*rec
603 END IF
604 END IF
605 x( j ) = zladiv( x( j ), tjjs )
606 xj = cabs1( x( j ) )
607 ELSE IF( tjj.GT.zero ) THEN
608*
609* 0 < abs(A(j,j)) <= SMLNUM:
610*
611 IF( xj.GT.tjj*bignum ) THEN
612*
613* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
614* to avoid overflow when dividing by A(j,j).
615*
616 rec = ( tjj*bignum ) / xj
617 IF( cnorm( j ).GT.one ) THEN
618*
619* Scale by 1/CNORM(j) to avoid overflow when
620* multiplying x(j) times column j.
621*
622 rec = rec / cnorm( j )
623 END IF
624 CALL zdscal( n, rec, x, 1 )
625 scale = scale*rec
626 xmax = xmax*rec
627 END IF
628 x( j ) = zladiv( x( j ), tjjs )
629 xj = cabs1( x( j ) )
630 ELSE
631*
632* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
633* scale = 0, and compute a solution to A*x = 0.
634*
635 DO 100 i = 1, n
636 x( i ) = zero
637 100 CONTINUE
638 x( j ) = one
639 xj = one
640 scale = zero
641 xmax = zero
642 END IF
643 110 CONTINUE
644*
645* Scale x if necessary to avoid overflow when adding a
646* multiple of column j of A.
647*
648 IF( xj.GT.one ) THEN
649 rec = one / xj
650 IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
651*
652* Scale x by 1/(2*abs(x(j))).
653*
654 rec = rec*half
655 CALL zdscal( n, rec, x, 1 )
656 scale = scale*rec
657 END IF
658 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
659*
660* Scale x by 1/2.
661*
662 CALL zdscal( n, half, x, 1 )
663 scale = scale*half
664 END IF
665*
666 IF( upper ) THEN
667 IF( j.GT.1 ) THEN
668*
669* Compute the update
670* x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
671* x(j)* A(max(1,j-kd):j-1,j)
672*
673 jlen = min( kd, j-1 )
674 CALL zaxpy( jlen, -x( j )*tscal,
675 \$ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
676 i = izamax( j-1, x, 1 )
677 xmax = cabs1( x( i ) )
678 END IF
679 ELSE IF( j.LT.n ) THEN
680*
681* Compute the update
682* x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
683* x(j) * A(j+1:min(j+kd,n),j)
684*
685 jlen = min( kd, n-j )
686 IF( jlen.GT.0 )
687 \$ CALL zaxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
688 \$ x( j+1 ), 1 )
689 i = j + izamax( n-j, x( j+1 ), 1 )
690 xmax = cabs1( x( i ) )
691 END IF
692 120 CONTINUE
693*
694 ELSE IF( lsame( trans, 'T' ) ) THEN
695*
696* Solve A**T * x = b
697*
698 DO 170 j = jfirst, jlast, jinc
699*
700* Compute x(j) = b(j) - sum A(k,j)*x(k).
701* k<>j
702*
703 xj = cabs1( x( j ) )
704 uscal = tscal
705 rec = one / max( xmax, one )
706 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
707*
708* If x(j) could overflow, scale x by 1/(2*XMAX).
709*
710 rec = rec*half
711 IF( nounit ) THEN
712 tjjs = ab( maind, j )*tscal
713 ELSE
714 tjjs = tscal
715 END IF
716 tjj = cabs1( tjjs )
717 IF( tjj.GT.one ) THEN
718*
719* Divide by A(j,j) when scaling x if A(j,j) > 1.
720*
721 rec = min( one, rec*tjj )
722 uscal = zladiv( uscal, tjjs )
723 END IF
724 IF( rec.LT.one ) THEN
725 CALL zdscal( n, rec, x, 1 )
726 scale = scale*rec
727 xmax = xmax*rec
728 END IF
729 END IF
730*
731 csumj = zero
732 IF( uscal.EQ.dcmplx( one ) ) THEN
733*
734* If the scaling needed for A in the dot product is 1,
735* call ZDOTU to perform the dot product.
736*
737 IF( upper ) THEN
738 jlen = min( kd, j-1 )
739 csumj = zdotu( jlen, ab( kd+1-jlen, j ), 1,
740 \$ x( j-jlen ), 1 )
741 ELSE
742 jlen = min( kd, n-j )
743 IF( jlen.GT.1 )
744 \$ csumj = zdotu( jlen, ab( 2, j ), 1, x( j+1 ),
745 \$ 1 )
746 END IF
747 ELSE
748*
749* Otherwise, use in-line code for the dot product.
750*
751 IF( upper ) THEN
752 jlen = min( kd, j-1 )
753 DO 130 i = 1, jlen
754 csumj = csumj + ( ab( kd+i-jlen, j )*uscal )*
755 \$ x( j-jlen-1+i )
756 130 CONTINUE
757 ELSE
758 jlen = min( kd, n-j )
759 DO 140 i = 1, jlen
760 csumj = csumj + ( ab( i+1, j )*uscal )*x( j+i )
761 140 CONTINUE
762 END IF
763 END IF
764*
765 IF( uscal.EQ.dcmplx( tscal ) ) THEN
766*
767* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
768* was not used to scale the dotproduct.
769*
770 x( j ) = x( j ) - csumj
771 xj = cabs1( x( j ) )
772 IF( nounit ) THEN
773*
774* Compute x(j) = x(j) / A(j,j), scaling if necessary.
775*
776 tjjs = ab( maind, j )*tscal
777 ELSE
778 tjjs = tscal
779 IF( tscal.EQ.one )
780 \$ GO TO 160
781 END IF
782 tjj = cabs1( tjjs )
783 IF( tjj.GT.smlnum ) THEN
784*
785* abs(A(j,j)) > SMLNUM:
786*
787 IF( tjj.LT.one ) THEN
788 IF( xj.GT.tjj*bignum ) THEN
789*
790* Scale X by 1/abs(x(j)).
791*
792 rec = one / xj
793 CALL zdscal( n, rec, x, 1 )
794 scale = scale*rec
795 xmax = xmax*rec
796 END IF
797 END IF
798 x( j ) = zladiv( x( j ), tjjs )
799 ELSE IF( tjj.GT.zero ) THEN
800*
801* 0 < abs(A(j,j)) <= SMLNUM:
802*
803 IF( xj.GT.tjj*bignum ) THEN
804*
805* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
806*
807 rec = ( tjj*bignum ) / xj
808 CALL zdscal( n, rec, x, 1 )
809 scale = scale*rec
810 xmax = xmax*rec
811 END IF
812 x( j ) = zladiv( x( j ), tjjs )
813 ELSE
814*
815* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
816* scale = 0 and compute a solution to A**T *x = 0.
817*
818 DO 150 i = 1, n
819 x( i ) = zero
820 150 CONTINUE
821 x( j ) = one
822 scale = zero
823 xmax = zero
824 END IF
825 160 CONTINUE
826 ELSE
827*
828* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
829* product has already been divided by 1/A(j,j).
830*
831 x( j ) = zladiv( x( j ), tjjs ) - csumj
832 END IF
833 xmax = max( xmax, cabs1( x( j ) ) )
834 170 CONTINUE
835*
836 ELSE
837*
838* Solve A**H * x = b
839*
840 DO 220 j = jfirst, jlast, jinc
841*
842* Compute x(j) = b(j) - sum A(k,j)*x(k).
843* k<>j
844*
845 xj = cabs1( x( j ) )
846 uscal = tscal
847 rec = one / max( xmax, one )
848 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
849*
850* If x(j) could overflow, scale x by 1/(2*XMAX).
851*
852 rec = rec*half
853 IF( nounit ) THEN
854 tjjs = dconjg( ab( maind, j ) )*tscal
855 ELSE
856 tjjs = tscal
857 END IF
858 tjj = cabs1( tjjs )
859 IF( tjj.GT.one ) THEN
860*
861* Divide by A(j,j) when scaling x if A(j,j) > 1.
862*
863 rec = min( one, rec*tjj )
864 uscal = zladiv( uscal, tjjs )
865 END IF
866 IF( rec.LT.one ) THEN
867 CALL zdscal( n, rec, x, 1 )
868 scale = scale*rec
869 xmax = xmax*rec
870 END IF
871 END IF
872*
873 csumj = zero
874 IF( uscal.EQ.dcmplx( one ) ) THEN
875*
876* If the scaling needed for A in the dot product is 1,
877* call ZDOTC to perform the dot product.
878*
879 IF( upper ) THEN
880 jlen = min( kd, j-1 )
881 csumj = zdotc( jlen, ab( kd+1-jlen, j ), 1,
882 \$ x( j-jlen ), 1 )
883 ELSE
884 jlen = min( kd, n-j )
885 IF( jlen.GT.1 )
886 \$ csumj = zdotc( jlen, ab( 2, j ), 1, x( j+1 ),
887 \$ 1 )
888 END IF
889 ELSE
890*
891* Otherwise, use in-line code for the dot product.
892*
893 IF( upper ) THEN
894 jlen = min( kd, j-1 )
895 DO 180 i = 1, jlen
896 csumj = csumj + ( dconjg( ab( kd+i-jlen, j ) )*
897 \$ uscal )*x( j-jlen-1+i )
898 180 CONTINUE
899 ELSE
900 jlen = min( kd, n-j )
901 DO 190 i = 1, jlen
902 csumj = csumj + ( dconjg( ab( i+1, j ) )*uscal )
903 \$ *x( j+i )
904 190 CONTINUE
905 END IF
906 END IF
907*
908 IF( uscal.EQ.dcmplx( tscal ) ) THEN
909*
910* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
911* was not used to scale the dotproduct.
912*
913 x( j ) = x( j ) - csumj
914 xj = cabs1( x( j ) )
915 IF( nounit ) THEN
916*
917* Compute x(j) = x(j) / A(j,j), scaling if necessary.
918*
919 tjjs = dconjg( ab( maind, j ) )*tscal
920 ELSE
921 tjjs = tscal
922 IF( tscal.EQ.one )
923 \$ GO TO 210
924 END IF
925 tjj = cabs1( tjjs )
926 IF( tjj.GT.smlnum ) THEN
927*
928* abs(A(j,j)) > SMLNUM:
929*
930 IF( tjj.LT.one ) THEN
931 IF( xj.GT.tjj*bignum ) THEN
932*
933* Scale X by 1/abs(x(j)).
934*
935 rec = one / xj
936 CALL zdscal( n, rec, x, 1 )
937 scale = scale*rec
938 xmax = xmax*rec
939 END IF
940 END IF
941 x( j ) = zladiv( x( j ), tjjs )
942 ELSE IF( tjj.GT.zero ) THEN
943*
944* 0 < abs(A(j,j)) <= SMLNUM:
945*
946 IF( xj.GT.tjj*bignum ) THEN
947*
948* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
949*
950 rec = ( tjj*bignum ) / xj
951 CALL zdscal( n, rec, x, 1 )
952 scale = scale*rec
953 xmax = xmax*rec
954 END IF
955 x( j ) = zladiv( x( j ), tjjs )
956 ELSE
957*
958* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
959* scale = 0 and compute a solution to A**H *x = 0.
960*
961 DO 200 i = 1, n
962 x( i ) = zero
963 200 CONTINUE
964 x( j ) = one
965 scale = zero
966 xmax = zero
967 END IF
968 210 CONTINUE
969 ELSE
970*
971* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
972* product has already been divided by 1/A(j,j).
973*
974 x( j ) = zladiv( x( j ), tjjs ) - csumj
975 END IF
976 xmax = max( xmax, cabs1( x( j ) ) )
977 220 CONTINUE
978 END IF
979 scale = scale / tscal
980 END IF
981*
982* Scale the column norms by 1/TSCAL for return.
983*
984 IF( tscal.NE.one ) THEN
985 CALL dscal( n, one / tscal, cnorm, 1 )
986 END IF
987*
988 RETURN
989*
990* End of ZLATBS
991*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
double precision function dzasum(n, zx, incx)
DZASUM
Definition dzasum.f:72
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
complex *16 function zdotc(n, zx, incx, zy, incy)
ZDOTC
Definition zdotc.f:83
complex *16 function zdotu(n, zx, incx, zy, incy)
ZDOTU
Definition zdotu.f:83
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
subroutine ztbsv(uplo, trans, diag, n, k, a, lda, x, incx)
ZTBSV
Definition ztbsv.f:189
Here is the call graph for this function:
Here is the caller graph for this function: