LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slatbs()

subroutine slatbs ( character  UPLO,
character  TRANS,
character  DIAG,
character  NORMIN,
integer  N,
integer  KD,
real, dimension( ldab, * )  AB,
integer  LDAB,
real, dimension( * )  X,
real  SCALE,
real, dimension( * )  CNORM,
integer  INFO 
)

SLATBS solves a triangular banded system of equations.

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

Purpose:
 SLATBS 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 STBSV 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 REAL 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 REAL 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  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 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, STBSV
  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 STBSV 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 STBSV if 1/M(n) and 1/G(n) are both greater
  than max(underflow, 1/overflow).

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