LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ clatbs()

subroutine clatbs ( character uplo,
character trans,
character diag,
character normin,
integer n,
integer kd,
complex, dimension( ldab, * ) ab,
integer ldab,
complex, dimension( * ) x,
real scale,
real, dimension( * ) cnorm,
integer info )

CLATBS solves a triangular banded system of equations.

Download CLATBS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> CLATBS 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 CTBSV 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 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 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 REAL
!>          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 REAL 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
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  A rough bound on x is computed; if that is less than overflow, CTBSV
!>  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 CTBSV 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 CTBSV if 1/M(n) and 1/G(n) are both greater
!>  than max(underflow, 1/overflow).
!>

Definition at line 239 of file clatbs.f.

242*
243* -- LAPACK auxiliary routine --
244* -- LAPACK is a software package provided by Univ. of Tennessee, --
245* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
246*
247* .. Scalar Arguments ..
248 CHARACTER DIAG, NORMIN, TRANS, UPLO
249 INTEGER INFO, KD, LDAB, N
250 REAL SCALE
251* ..
252* .. Array Arguments ..
253 REAL CNORM( * )
254 COMPLEX AB( LDAB, * ), X( * )
255* ..
256*
257* =====================================================================
258*
259* .. Parameters ..
260 REAL ZERO, HALF, ONE, TWO
261 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
262 $ two = 2.0e+0 )
263* ..
264* .. Local Scalars ..
265 LOGICAL NOTRAN, NOUNIT, UPPER
266 INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
267 REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
268 $ XBND, XJ, XMAX
269 COMPLEX CSUMJ, TJJS, USCAL, ZDUM
270* ..
271* .. External Functions ..
272 LOGICAL LSAME
273 INTEGER ICAMAX, ISAMAX
274 REAL SCASUM, SLAMCH
275 COMPLEX CDOTC, CDOTU, CLADIV
276 EXTERNAL lsame, icamax, isamax, scasum, slamch,
277 $ cdotc,
278 $ cdotu, cladiv
279* ..
280* .. External Subroutines ..
281 EXTERNAL caxpy, csscal, ctbsv, sscal, xerbla
282* ..
283* .. Intrinsic Functions ..
284 INTRINSIC abs, aimag, cmplx, conjg, max, min, real
285* ..
286* .. Statement Functions ..
287 REAL CABS1, CABS2
288* ..
289* .. Statement Function definitions ..
290 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
291 cabs2( zdum ) = abs( real( zdum ) / 2. ) +
292 $ abs( aimag( zdum ) / 2. )
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( 'CLATBS', -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 = slamch( 'Safe minimum' ) / slamch( '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 ) = scasum( 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 ) = scasum( 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 = isamax( 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 sscal( 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 CTBSV 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 ctbsv( 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 csscal( 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 110 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 105
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 csscal( n, rec, x, 1 )
601 scale = scale*rec
602 xmax = xmax*rec
603 END IF
604 END IF
605 x( j ) = cladiv( 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 csscal( n, rec, x, 1 )
625 scale = scale*rec
626 xmax = xmax*rec
627 END IF
628 x( j ) = cladiv( 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 105 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 csscal( 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 csscal( 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 caxpy( jlen, -x( j )*tscal,
675 $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
676 i = icamax( 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 caxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
688 $ x( j+1 ), 1 )
689 i = j + icamax( n-j, x( j+1 ), 1 )
690 xmax = cabs1( x( i ) )
691 END IF
692 110 CONTINUE
693*
694 ELSE IF( lsame( trans, 'T' ) ) THEN
695*
696* Solve A**T * x = b
697*
698 DO 150 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 = cladiv( uscal, tjjs )
723 END IF
724 IF( rec.LT.one ) THEN
725 CALL csscal( 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.cmplx( one ) ) THEN
733*
734* If the scaling needed for A in the dot product is 1,
735* call CDOTU to perform the dot product.
736*
737 IF( upper ) THEN
738 jlen = min( kd, j-1 )
739 csumj = cdotu( 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 = cdotu( 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 120 i = 1, jlen
754 csumj = csumj + ( ab( kd+i-jlen, j )*uscal )*
755 $ x( j-jlen-1+i )
756 120 CONTINUE
757 ELSE
758 jlen = min( kd, n-j )
759 DO 130 i = 1, jlen
760 csumj = csumj + ( ab( i+1, j )*uscal )*x( j+i )
761 130 CONTINUE
762 END IF
763 END IF
764*
765 IF( uscal.EQ.cmplx( 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 145
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 csscal( n, rec, x, 1 )
794 scale = scale*rec
795 xmax = xmax*rec
796 END IF
797 END IF
798 x( j ) = cladiv( 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 csscal( n, rec, x, 1 )
809 scale = scale*rec
810 xmax = xmax*rec
811 END IF
812 x( j ) = cladiv( 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 140 i = 1, n
819 x( i ) = zero
820 140 CONTINUE
821 x( j ) = one
822 scale = zero
823 xmax = zero
824 END IF
825 145 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 ) = cladiv( x( j ), tjjs ) - csumj
832 END IF
833 xmax = max( xmax, cabs1( x( j ) ) )
834 150 CONTINUE
835*
836 ELSE
837*
838* Solve A**H * x = b
839*
840 DO 190 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 = conjg( 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 = cladiv( uscal, tjjs )
865 END IF
866 IF( rec.LT.one ) THEN
867 CALL csscal( 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.cmplx( one ) ) THEN
875*
876* If the scaling needed for A in the dot product is 1,
877* call CDOTC to perform the dot product.
878*
879 IF( upper ) THEN
880 jlen = min( kd, j-1 )
881 csumj = cdotc( 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 = cdotc( 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 160 i = 1, jlen
896 csumj = csumj + ( conjg( ab( kd+i-jlen, j ) )*
897 $ uscal )*x( j-jlen-1+i )
898 160 CONTINUE
899 ELSE
900 jlen = min( kd, n-j )
901 DO 170 i = 1, jlen
902 csumj = csumj + ( conjg( ab( i+1, j ) )*uscal )*
903 $ x( j+i )
904 170 CONTINUE
905 END IF
906 END IF
907*
908 IF( uscal.EQ.cmplx( 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 = conjg( ab( maind, j ) )*tscal
920 ELSE
921 tjjs = tscal
922 IF( tscal.EQ.one )
923 $ GO TO 185
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 csscal( n, rec, x, 1 )
937 scale = scale*rec
938 xmax = xmax*rec
939 END IF
940 END IF
941 x( j ) = cladiv( 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 csscal( n, rec, x, 1 )
952 scale = scale*rec
953 xmax = xmax*rec
954 END IF
955 x( j ) = cladiv( 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 180 i = 1, n
962 x( i ) = zero
963 180 CONTINUE
964 x( j ) = one
965 scale = zero
966 xmax = zero
967 END IF
968 185 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 ) = cladiv( x( j ), tjjs ) - csumj
975 END IF
976 xmax = max( xmax, cabs1( x( j ) ) )
977 190 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 sscal( n, one / tscal, cnorm, 1 )
986 END IF
987*
988 RETURN
989*
990* End of CLATBS
991*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scasum
real function scasum(n, cx, incx)
SCASUM
Definition scasum.f:72
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
cdotu
complex function cdotu(n, cx, incx, cy, incy)
CDOTU
Definition cdotu.f:83
cdotc
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
icamax
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
cladiv
complex function cladiv(x, y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition cladiv.f:62
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
ctbsv
subroutine ctbsv(uplo, trans, diag, n, k, a, lda, x, incx)
CTBSV
Definition ctbsv.f:189
Here is the call graph for this function:
Here is the caller graph for this function: