LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
slatbs.f
Go to the documentation of this file.
1 *> \brief \b SLATBS solves a triangular banded system of equations.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SLATBS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slatbs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slatbs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slatbs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
22 * SCALE, CNORM, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER DIAG, NORMIN, TRANS, UPLO
26 * INTEGER INFO, KD, LDAB, N
27 * REAL SCALE
28 * ..
29 * .. Array Arguments ..
30 * REAL AB( LDAB, * ), CNORM( * ), X( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SLATBS solves one of the triangular systems
40 *>
41 *> A *x = s*b or A**T*x = s*b
42 *>
43 *> with scaling to prevent overflow, where A is an upper or lower
44 *> triangular band matrix. Here A**T denotes the transpose of A, x and b
45 *> are n-element vectors, and s is a scaling factor, usually less than
46 *> or equal to 1, chosen so that the components of x will be less than
47 *> the overflow threshold. If the unscaled problem will not cause
48 *> overflow, the Level 2 BLAS routine STBSV is called. If the matrix A
49 *> is singular (A(j,j) = 0 for some j), then s is set to 0 and a
50 *> non-trivial solution to A*x = 0 is returned.
51 *> \endverbatim
52 *
53 * Arguments:
54 * ==========
55 *
56 *> \param[in] UPLO
57 *> \verbatim
58 *> UPLO is CHARACTER*1
59 *> Specifies whether the matrix A is upper or lower triangular.
60 *> = 'U': Upper triangular
61 *> = 'L': Lower triangular
62 *> \endverbatim
63 *>
64 *> \param[in] TRANS
65 *> \verbatim
66 *> TRANS is CHARACTER*1
67 *> Specifies the operation applied to A.
68 *> = 'N': Solve A * x = s*b (No transpose)
69 *> = 'T': Solve A**T* x = s*b (Transpose)
70 *> = 'C': Solve A**T* x = s*b (Conjugate transpose = Transpose)
71 *> \endverbatim
72 *>
73 *> \param[in] DIAG
74 *> \verbatim
75 *> DIAG is CHARACTER*1
76 *> Specifies whether or not the matrix A is unit triangular.
77 *> = 'N': Non-unit triangular
78 *> = 'U': Unit triangular
79 *> \endverbatim
80 *>
81 *> \param[in] NORMIN
82 *> \verbatim
83 *> NORMIN is CHARACTER*1
84 *> Specifies whether CNORM has been set or not.
85 *> = 'Y': CNORM contains the column norms on entry
86 *> = 'N': CNORM is not set on entry. On exit, the norms will
87 *> be computed and stored in CNORM.
88 *> \endverbatim
89 *>
90 *> \param[in] N
91 *> \verbatim
92 *> N is INTEGER
93 *> The order of the matrix A. N >= 0.
94 *> \endverbatim
95 *>
96 *> \param[in] KD
97 *> \verbatim
98 *> KD is INTEGER
99 *> The number of subdiagonals or superdiagonals in the
100 *> triangular matrix A. KD >= 0.
101 *> \endverbatim
102 *>
103 *> \param[in] AB
104 *> \verbatim
105 *> AB is REAL array, dimension (LDAB,N)
106 *> The upper or lower triangular band matrix A, stored in the
107 *> first KD+1 rows of the array. The j-th column of A is stored
108 *> in the j-th column of the array AB as follows:
109 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
110 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
111 *> \endverbatim
112 *>
113 *> \param[in] LDAB
114 *> \verbatim
115 *> LDAB is INTEGER
116 *> The leading dimension of the array AB. LDAB >= KD+1.
117 *> \endverbatim
118 *>
119 *> \param[in,out] X
120 *> \verbatim
121 *> X is REAL array, dimension (N)
122 *> On entry, the right hand side b of the triangular system.
123 *> On exit, X is overwritten by the solution vector x.
124 *> \endverbatim
125 *>
126 *> \param[out] SCALE
127 *> \verbatim
128 *> SCALE is REAL
129 *> The scaling factor s for the triangular system
130 *> A * x = s*b or A**T* x = s*b.
131 *> If SCALE = 0, the matrix A is singular or badly scaled, and
132 *> the vector x is an exact or approximate solution to A*x = 0.
133 *> \endverbatim
134 *>
135 *> \param[in,out] CNORM
136 *> \verbatim
137 *> CNORM is REAL array, dimension (N)
138 *>
139 *> If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
140 *> contains the norm of the off-diagonal part of the j-th column
141 *> of A. If TRANS = 'N', CNORM(j) must be greater than or equal
142 *> to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
143 *> must be greater than or equal to the 1-norm.
144 *>
145 *> If NORMIN = 'N', CNORM is an output argument and CNORM(j)
146 *> returns the 1-norm of the offdiagonal part of the j-th column
147 *> of A.
148 *> \endverbatim
149 *>
150 *> \param[out] INFO
151 *> \verbatim
152 *> INFO is INTEGER
153 *> = 0: successful exit
154 *> < 0: if INFO = -k, the k-th argument had an illegal value
155 *> \endverbatim
156 *
157 * Authors:
158 * ========
159 *
160 *> \author Univ. of Tennessee
161 *> \author Univ. of California Berkeley
162 *> \author Univ. of Colorado Denver
163 *> \author NAG Ltd.
164 *
165 *> \date September 2012
166 *
167 *> \ingroup realOTHERauxiliary
168 *
169 *> \par Further Details:
170 * =====================
171 *>
172 *> \verbatim
173 *>
174 *> A rough bound on x is computed; if that is less than overflow, STBSV
175 *> is called, otherwise, specific code is used which checks for possible
176 *> overflow or divide-by-zero at every operation.
177 *>
178 *> A columnwise scheme is used for solving A*x = b. The basic algorithm
179 *> if A is lower triangular is
180 *>
181 *> x[1:n] := b[1:n]
182 *> for j = 1, ..., n
183 *> x(j) := x(j) / A(j,j)
184 *> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
185 *> end
186 *>
187 *> Define bounds on the components of x after j iterations of the loop:
188 *> M(j) = bound on x[1:j]
189 *> G(j) = bound on x[j+1:n]
190 *> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
191 *>
192 *> Then for iteration j+1 we have
193 *> M(j+1) <= G(j) / | A(j+1,j+1) |
194 *> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
195 *> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
196 *>
197 *> where CNORM(j+1) is greater than or equal to the infinity-norm of
198 *> column j+1 of A, not counting the diagonal. Hence
199 *>
200 *> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
201 *> 1<=i<=j
202 *> and
203 *>
204 *> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
205 *> 1<=i< j
206 *>
207 *> Since |x(j)| <= M(j), we use the Level 2 BLAS routine STBSV if the
208 *> reciprocal of the largest M(j), j=1,..,n, is larger than
209 *> max(underflow, 1/overflow).
210 *>
211 *> The bound on x(j) is also used to determine when a step in the
212 *> columnwise method can be performed without fear of overflow. If
213 *> the computed bound is greater than a large constant, x is scaled to
214 *> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
215 *> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
216 *>
217 *> Similarly, a row-wise scheme is used to solve A**T*x = b. The basic
218 *> algorithm for A upper triangular is
219 *>
220 *> for j = 1, ..., n
221 *> x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
222 *> end
223 *>
224 *> We simultaneously compute two bounds
225 *> G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
226 *> M(j) = bound on x(i), 1<=i<=j
227 *>
228 *> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
229 *> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
230 *> Then the bound on x(j) is
231 *>
232 *> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
233 *>
234 *> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
235 *> 1<=i<=j
236 *>
237 *> and we can safely call STBSV if 1/M(n) and 1/G(n) are both greater
238 *> than max(underflow, 1/overflow).
239 *> \endverbatim
240 *>
241 * =====================================================================
242  SUBROUTINE slatbs( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
243  $ scale, cnorm, info )
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  REAL SCALE
254 * ..
255 * .. Array Arguments ..
256  REAL AB( ldab, * ), CNORM( * ), X( * )
257 * ..
258 *
259 * =====================================================================
260 *
261 * .. Parameters ..
262  REAL ZERO, HALF, ONE
263  parameter ( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0 )
264 * ..
265 * .. Local Scalars ..
266  LOGICAL NOTRAN, NOUNIT, UPPER
267  INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
268  REAL BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
269  $ tmax, tscal, uscal, xbnd, xj, xmax
270 * ..
271 * .. External Functions ..
272  LOGICAL LSAME
273  INTEGER ISAMAX
274  REAL SASUM, SDOT, SLAMCH
275  EXTERNAL lsame, isamax, sasum, sdot, slamch
276 * ..
277 * .. External Subroutines ..
278  EXTERNAL saxpy, sscal, stbsv, 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( 'SLATBS', -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 = slamch( 'Safe minimum' ) / slamch( '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 ) = sasum( 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 ) = sasum( 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 = isamax( 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 sscal( 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 STBSV can be used.
366 *
367  j = isamax( 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 stbsv( 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 sscal( n, scale, x, 1 )
534  xmax = bignum
535  END IF
536 *
537  IF( notran ) THEN
538 *
539 * Solve A * x = b
540 *
541  DO 100 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 95
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 sscal( 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 sscal( 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  95 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 sscal( 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 sscal( 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 saxpy( jlen, -x( j )*tscal,
639  $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
640  i = isamax( 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 saxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
652  $ x( j+1 ), 1 )
653  i = j + isamax( n-j, x( j+1 ), 1 )
654  xmax = abs( x( i ) )
655  END IF
656  100 CONTINUE
657 *
658  ELSE
659 *
660 * Solve A**T * x = b
661 *
662  DO 140 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 sscal( 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 SDOT to perform the dot product.
700 *
701  IF( upper ) THEN
702  jlen = min( kd, j-1 )
703  sumj = sdot( 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 = sdot( 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 110 i = 1, jlen
717  sumj = sumj + ( ab( kd+i-jlen, j )*uscal )*
718  $ x( j-jlen-1+i )
719  110 CONTINUE
720  ELSE
721  jlen = min( kd, n-j )
722  DO 120 i = 1, jlen
723  sumj = sumj + ( ab( i+1, j )*uscal )*x( j+i )
724  120 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 135
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 sscal( 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 sscal( 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 130 i = 1, n
782  x( i ) = zero
783  130 CONTINUE
784  x( j ) = one
785  scale = zero
786  xmax = zero
787  END IF
788  135 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  140 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 sscal( n, one / tscal, cnorm, 1 )
806  END IF
807 *
808  RETURN
809 *
810 * End of SLATBS
811 *
812  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
stbsv
subroutine stbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
STBSV
Definition: stbsv.f:191
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
slatbs
subroutine slatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
SLATBS solves a triangular banded system of equations.
Definition: slatbs.f:244