LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlatbs ( character  UPLO,
character  TRANS,
character  DIAG,
character  NORMIN,
integer  N,
integer  KD,
double precision, dimension( ldab, * )  AB,
integer  LDAB,
double precision, dimension( * )  X,
double precision  SCALE,
double precision, dimension( * )  CNORM,
integer  INFO 
)

DLATBS solves a triangular banded system of equations.

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

Purpose:
 DLATBS solves one of the triangular systems

    A *x = s*b  or  A**T*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 DTBSV 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**T* x = s*b  (Conjugate transpose = 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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  or  A**T* 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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  A rough bound on x is computed; if that is less than overflow, DTBSV
  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 DTBSV 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.  The basic
  algorithm for A upper triangular is

       for j = 1, ..., n
            x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
       end

  We simultaneously compute two bounds
       G(j) = bound on ( b(i) - A[1:i-1,i]**T * 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 DTBSV if 1/M(n) and 1/G(n) are both greater
  than max(underflow, 1/overflow).

Definition at line 244 of file dlatbs.f.

244 *
245 * -- LAPACK auxiliary routine (version 3.4.2) --
246 * -- LAPACK is a software package provided by Univ. of Tennessee, --
247 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
248 * September 2012
249 *
250 * .. Scalar Arguments ..
251  CHARACTER diag, normin, trans, uplo
252  INTEGER info, kd, ldab, n
253  DOUBLE PRECISION scale
254 * ..
255 * .. Array Arguments ..
256  DOUBLE PRECISION ab( ldab, * ), cnorm( * ), x( * )
257 * ..
258 *
259 * =====================================================================
260 *
261 * .. Parameters ..
262  DOUBLE PRECISION zero, half, one
263  parameter ( zero = 0.0d+0, half = 0.5d+0, one = 1.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, sumj, tjj, tjjs,
269  $ tmax, tscal, uscal, xbnd, xj, xmax
270 * ..
271 * .. External Functions ..
272  LOGICAL lsame
273  INTEGER idamax
274  DOUBLE PRECISION dasum, ddot, dlamch
275  EXTERNAL lsame, idamax, dasum, ddot, dlamch
276 * ..
277 * .. External Subroutines ..
278  EXTERNAL daxpy, dscal, dtbsv, xerbla
279 * ..
280 * .. Intrinsic Functions ..
281  INTRINSIC abs, max, min
282 * ..
283 * .. Executable Statements ..
284 *
285  info = 0
286  upper = lsame( uplo, 'U' )
287  notran = lsame( trans, 'N' )
288  nounit = lsame( diag, 'N' )
289 *
290 * Test the input parameters.
291 *
292  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
293  info = -1
294  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
295  $ lsame( trans, 'C' ) ) THEN
296  info = -2
297  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
298  info = -3
299  ELSE IF( .NOT.lsame( normin, 'Y' ) .AND. .NOT.
300  $ lsame( normin, 'N' ) ) THEN
301  info = -4
302  ELSE IF( n.LT.0 ) THEN
303  info = -5
304  ELSE IF( kd.LT.0 ) THEN
305  info = -6
306  ELSE IF( ldab.LT.kd+1 ) THEN
307  info = -8
308  END IF
309  IF( info.NE.0 ) THEN
310  CALL xerbla( 'DLATBS', -info )
311  RETURN
312  END IF
313 *
314 * Quick return if possible
315 *
316  IF( n.EQ.0 )
317  $ RETURN
318 *
319 * Determine machine dependent parameters to control overflow.
320 *
321  smlnum = dlamch( 'Safe minimum' ) / dlamch( 'Precision' )
322  bignum = one / smlnum
323  scale = one
324 *
325  IF( lsame( normin, 'N' ) ) THEN
326 *
327 * Compute the 1-norm of each column, not including the diagonal.
328 *
329  IF( upper ) THEN
330 *
331 * A is upper triangular.
332 *
333  DO 10 j = 1, n
334  jlen = min( kd, j-1 )
335  cnorm( j ) = dasum( jlen, ab( kd+1-jlen, j ), 1 )
336  10 CONTINUE
337  ELSE
338 *
339 * A is lower triangular.
340 *
341  DO 20 j = 1, n
342  jlen = min( kd, n-j )
343  IF( jlen.GT.0 ) THEN
344  cnorm( j ) = dasum( jlen, ab( 2, j ), 1 )
345  ELSE
346  cnorm( j ) = zero
347  END IF
348  20 CONTINUE
349  END IF
350  END IF
351 *
352 * Scale the column norms by TSCAL if the maximum element in CNORM is
353 * greater than BIGNUM.
354 *
355  imax = idamax( n, cnorm, 1 )
356  tmax = cnorm( imax )
357  IF( tmax.LE.bignum ) THEN
358  tscal = one
359  ELSE
360  tscal = one / ( smlnum*tmax )
361  CALL dscal( n, tscal, cnorm, 1 )
362  END IF
363 *
364 * Compute a bound on the computed solution vector to see if the
365 * Level 2 BLAS routine DTBSV can be used.
366 *
367  j = idamax( n, x, 1 )
368  xmax = abs( x( j ) )
369  xbnd = xmax
370  IF( notran ) THEN
371 *
372 * Compute the growth in A * x = b.
373 *
374  IF( upper ) THEN
375  jfirst = n
376  jlast = 1
377  jinc = -1
378  maind = kd + 1
379  ELSE
380  jfirst = 1
381  jlast = n
382  jinc = 1
383  maind = 1
384  END IF
385 *
386  IF( tscal.NE.one ) THEN
387  grow = zero
388  GO TO 50
389  END IF
390 *
391  IF( nounit ) THEN
392 *
393 * A is non-unit triangular.
394 *
395 * Compute GROW = 1/G(j) and XBND = 1/M(j).
396 * Initially, G(0) = max{x(i), i=1,...,n}.
397 *
398  grow = one / max( xbnd, smlnum )
399  xbnd = grow
400  DO 30 j = jfirst, jlast, jinc
401 *
402 * Exit the loop if the growth factor is too small.
403 *
404  IF( grow.LE.smlnum )
405  $ GO TO 50
406 *
407 * M(j) = G(j-1) / abs(A(j,j))
408 *
409  tjj = abs( ab( maind, j ) )
410  xbnd = min( xbnd, min( one, tjj )*grow )
411  IF( tjj+cnorm( j ).GE.smlnum ) THEN
412 *
413 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
414 *
415  grow = grow*( tjj / ( tjj+cnorm( j ) ) )
416  ELSE
417 *
418 * G(j) could overflow, set GROW to 0.
419 *
420  grow = zero
421  END IF
422  30 CONTINUE
423  grow = xbnd
424  ELSE
425 *
426 * A is unit triangular.
427 *
428 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
429 *
430  grow = min( one, one / max( xbnd, smlnum ) )
431  DO 40 j = jfirst, jlast, jinc
432 *
433 * Exit the loop if the growth factor is too small.
434 *
435  IF( grow.LE.smlnum )
436  $ GO TO 50
437 *
438 * G(j) = G(j-1)*( 1 + CNORM(j) )
439 *
440  grow = grow*( one / ( one+cnorm( j ) ) )
441  40 CONTINUE
442  END IF
443  50 CONTINUE
444 *
445  ELSE
446 *
447 * Compute the growth in A**T * x = b.
448 *
449  IF( upper ) THEN
450  jfirst = 1
451  jlast = n
452  jinc = 1
453  maind = kd + 1
454  ELSE
455  jfirst = n
456  jlast = 1
457  jinc = -1
458  maind = 1
459  END IF
460 *
461  IF( tscal.NE.one ) THEN
462  grow = zero
463  GO TO 80
464  END IF
465 *
466  IF( nounit ) THEN
467 *
468 * A is non-unit triangular.
469 *
470 * Compute GROW = 1/G(j) and XBND = 1/M(j).
471 * Initially, M(0) = max{x(i), i=1,...,n}.
472 *
473  grow = one / max( xbnd, smlnum )
474  xbnd = grow
475  DO 60 j = jfirst, jlast, jinc
476 *
477 * Exit the loop if the growth factor is too small.
478 *
479  IF( grow.LE.smlnum )
480  $ GO TO 80
481 *
482 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
483 *
484  xj = one + cnorm( j )
485  grow = min( grow, xbnd / xj )
486 *
487 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
488 *
489  tjj = abs( ab( maind, j ) )
490  IF( xj.GT.tjj )
491  $ xbnd = xbnd*( tjj / xj )
492  60 CONTINUE
493  grow = min( grow, xbnd )
494  ELSE
495 *
496 * A is unit triangular.
497 *
498 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
499 *
500  grow = min( one, one / max( xbnd, smlnum ) )
501  DO 70 j = jfirst, jlast, jinc
502 *
503 * Exit the loop if the growth factor is too small.
504 *
505  IF( grow.LE.smlnum )
506  $ GO TO 80
507 *
508 * G(j) = ( 1 + CNORM(j) )*G(j-1)
509 *
510  xj = one + cnorm( j )
511  grow = grow / xj
512  70 CONTINUE
513  END IF
514  80 CONTINUE
515  END IF
516 *
517  IF( ( grow*tscal ).GT.smlnum ) THEN
518 *
519 * Use the Level 2 BLAS solve if the reciprocal of the bound on
520 * elements of X is not too small.
521 *
522  CALL dtbsv( uplo, trans, diag, n, kd, ab, ldab, x, 1 )
523  ELSE
524 *
525 * Use a Level 1 BLAS solve, scaling intermediate results.
526 *
527  IF( xmax.GT.bignum ) THEN
528 *
529 * Scale X so that its components are less than or equal to
530 * BIGNUM in absolute value.
531 *
532  scale = bignum / xmax
533  CALL dscal( n, scale, x, 1 )
534  xmax = bignum
535  END IF
536 *
537  IF( notran ) THEN
538 *
539 * Solve A * x = b
540 *
541  DO 110 j = jfirst, jlast, jinc
542 *
543 * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
544 *
545  xj = abs( x( j ) )
546  IF( nounit ) THEN
547  tjjs = ab( maind, j )*tscal
548  ELSE
549  tjjs = tscal
550  IF( tscal.EQ.one )
551  $ GO TO 100
552  END IF
553  tjj = abs( tjjs )
554  IF( tjj.GT.smlnum ) THEN
555 *
556 * abs(A(j,j)) > SMLNUM:
557 *
558  IF( tjj.LT.one ) THEN
559  IF( xj.GT.tjj*bignum ) THEN
560 *
561 * Scale x by 1/b(j).
562 *
563  rec = one / xj
564  CALL dscal( n, rec, x, 1 )
565  scale = scale*rec
566  xmax = xmax*rec
567  END IF
568  END IF
569  x( j ) = x( j ) / tjjs
570  xj = abs( x( j ) )
571  ELSE IF( tjj.GT.zero ) THEN
572 *
573 * 0 < abs(A(j,j)) <= SMLNUM:
574 *
575  IF( xj.GT.tjj*bignum ) THEN
576 *
577 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
578 * to avoid overflow when dividing by A(j,j).
579 *
580  rec = ( tjj*bignum ) / xj
581  IF( cnorm( j ).GT.one ) THEN
582 *
583 * Scale by 1/CNORM(j) to avoid overflow when
584 * multiplying x(j) times column j.
585 *
586  rec = rec / cnorm( j )
587  END IF
588  CALL dscal( n, rec, x, 1 )
589  scale = scale*rec
590  xmax = xmax*rec
591  END IF
592  x( j ) = x( j ) / tjjs
593  xj = abs( x( j ) )
594  ELSE
595 *
596 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
597 * scale = 0, and compute a solution to A*x = 0.
598 *
599  DO 90 i = 1, n
600  x( i ) = zero
601  90 CONTINUE
602  x( j ) = one
603  xj = one
604  scale = zero
605  xmax = zero
606  END IF
607  100 CONTINUE
608 *
609 * Scale x if necessary to avoid overflow when adding a
610 * multiple of column j of A.
611 *
612  IF( xj.GT.one ) THEN
613  rec = one / xj
614  IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
615 *
616 * Scale x by 1/(2*abs(x(j))).
617 *
618  rec = rec*half
619  CALL dscal( n, rec, x, 1 )
620  scale = scale*rec
621  END IF
622  ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
623 *
624 * Scale x by 1/2.
625 *
626  CALL dscal( n, half, x, 1 )
627  scale = scale*half
628  END IF
629 *
630  IF( upper ) THEN
631  IF( j.GT.1 ) THEN
632 *
633 * Compute the update
634 * x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
635 * x(j)* A(max(1,j-kd):j-1,j)
636 *
637  jlen = min( kd, j-1 )
638  CALL daxpy( jlen, -x( j )*tscal,
639  $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
640  i = idamax( j-1, x, 1 )
641  xmax = abs( x( i ) )
642  END IF
643  ELSE IF( j.LT.n ) THEN
644 *
645 * Compute the update
646 * x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
647 * x(j) * A(j+1:min(j+kd,n),j)
648 *
649  jlen = min( kd, n-j )
650  IF( jlen.GT.0 )
651  $ CALL daxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
652  $ x( j+1 ), 1 )
653  i = j + idamax( n-j, x( j+1 ), 1 )
654  xmax = abs( x( i ) )
655  END IF
656  110 CONTINUE
657 *
658  ELSE
659 *
660 * Solve A**T * x = b
661 *
662  DO 160 j = jfirst, jlast, jinc
663 *
664 * Compute x(j) = b(j) - sum A(k,j)*x(k).
665 * k<>j
666 *
667  xj = abs( x( j ) )
668  uscal = tscal
669  rec = one / max( xmax, one )
670  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
671 *
672 * If x(j) could overflow, scale x by 1/(2*XMAX).
673 *
674  rec = rec*half
675  IF( nounit ) THEN
676  tjjs = ab( maind, j )*tscal
677  ELSE
678  tjjs = tscal
679  END IF
680  tjj = abs( tjjs )
681  IF( tjj.GT.one ) THEN
682 *
683 * Divide by A(j,j) when scaling x if A(j,j) > 1.
684 *
685  rec = min( one, rec*tjj )
686  uscal = uscal / tjjs
687  END IF
688  IF( rec.LT.one ) THEN
689  CALL dscal( n, rec, x, 1 )
690  scale = scale*rec
691  xmax = xmax*rec
692  END IF
693  END IF
694 *
695  sumj = zero
696  IF( uscal.EQ.one ) THEN
697 *
698 * If the scaling needed for A in the dot product is 1,
699 * call DDOT to perform the dot product.
700 *
701  IF( upper ) THEN
702  jlen = min( kd, j-1 )
703  sumj = ddot( jlen, ab( kd+1-jlen, j ), 1,
704  $ x( j-jlen ), 1 )
705  ELSE
706  jlen = min( kd, n-j )
707  IF( jlen.GT.0 )
708  $ sumj = ddot( jlen, ab( 2, j ), 1, x( j+1 ), 1 )
709  END IF
710  ELSE
711 *
712 * Otherwise, use in-line code for the dot product.
713 *
714  IF( upper ) THEN
715  jlen = min( kd, j-1 )
716  DO 120 i = 1, jlen
717  sumj = sumj + ( ab( kd+i-jlen, j )*uscal )*
718  $ x( j-jlen-1+i )
719  120 CONTINUE
720  ELSE
721  jlen = min( kd, n-j )
722  DO 130 i = 1, jlen
723  sumj = sumj + ( ab( i+1, j )*uscal )*x( j+i )
724  130 CONTINUE
725  END IF
726  END IF
727 *
728  IF( uscal.EQ.tscal ) THEN
729 *
730 * Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
731 * was not used to scale the dotproduct.
732 *
733  x( j ) = x( j ) - sumj
734  xj = abs( x( j ) )
735  IF( nounit ) THEN
736 *
737 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
738 *
739  tjjs = ab( maind, j )*tscal
740  ELSE
741  tjjs = tscal
742  IF( tscal.EQ.one )
743  $ GO TO 150
744  END IF
745  tjj = abs( tjjs )
746  IF( tjj.GT.smlnum ) THEN
747 *
748 * abs(A(j,j)) > SMLNUM:
749 *
750  IF( tjj.LT.one ) THEN
751  IF( xj.GT.tjj*bignum ) THEN
752 *
753 * Scale X by 1/abs(x(j)).
754 *
755  rec = one / xj
756  CALL dscal( n, rec, x, 1 )
757  scale = scale*rec
758  xmax = xmax*rec
759  END IF
760  END IF
761  x( j ) = x( j ) / tjjs
762  ELSE IF( tjj.GT.zero ) THEN
763 *
764 * 0 < abs(A(j,j)) <= SMLNUM:
765 *
766  IF( xj.GT.tjj*bignum ) THEN
767 *
768 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
769 *
770  rec = ( tjj*bignum ) / xj
771  CALL dscal( n, rec, x, 1 )
772  scale = scale*rec
773  xmax = xmax*rec
774  END IF
775  x( j ) = x( j ) / tjjs
776  ELSE
777 *
778 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
779 * scale = 0, and compute a solution to A**T*x = 0.
780 *
781  DO 140 i = 1, n
782  x( i ) = zero
783  140 CONTINUE
784  x( j ) = one
785  scale = zero
786  xmax = zero
787  END IF
788  150 CONTINUE
789  ELSE
790 *
791 * Compute x(j) := x(j) / A(j,j) - sumj if the dot
792 * product has already been divided by 1/A(j,j).
793 *
794  x( j ) = x( j ) / tjjs - sumj
795  END IF
796  xmax = max( xmax, abs( x( j ) ) )
797  160 CONTINUE
798  END IF
799  scale = scale / tscal
800  END IF
801 *
802 * Scale the column norms by 1/TSCAL for return.
803 *
804  IF( tscal.NE.one ) THEN
805  CALL dscal( n, one / tscal, cnorm, 1 )
806  END IF
807 *
808  RETURN
809 *
810 * End of DLATBS
811 *
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
subroutine dtbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
DTBSV
Definition: dtbsv.f:191
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
double precision function dasum(N, DX, INCX)
DASUM
Definition: dasum.f:53
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: