LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlatbs.f
Go to the documentation of this file.
1 *> \brief \b DLATBS 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 DLATBS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlatbs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlatbs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlatbs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLATBS( 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 * DOUBLE PRECISION SCALE
28 * ..
29 * .. Array Arguments ..
30 * DOUBLE PRECISION AB( LDAB, * ), CNORM( * ), X( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DLATBS 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 DTBSV 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 *> \ingroup doubleOTHERauxiliary
166 *
167 *> \par Further Details:
168 * =====================
169 *>
170 *> \verbatim
171 *>
172 *> A rough bound on x is computed; if that is less than overflow, DTBSV
173 *> is called, otherwise, specific code is used which checks for possible
174 *> overflow or divide-by-zero at every operation.
175 *>
176 *> A columnwise scheme is used for solving A*x = b. The basic algorithm
177 *> if A is lower triangular is
178 *>
179 *> x[1:n] := b[1:n]
180 *> for j = 1, ..., n
181 *> x(j) := x(j) / A(j,j)
182 *> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
183 *> end
184 *>
185 *> Define bounds on the components of x after j iterations of the loop:
186 *> M(j) = bound on x[1:j]
187 *> G(j) = bound on x[j+1:n]
188 *> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
189 *>
190 *> Then for iteration j+1 we have
191 *> M(j+1) <= G(j) / | A(j+1,j+1) |
192 *> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
193 *> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
194 *>
195 *> where CNORM(j+1) is greater than or equal to the infinity-norm of
196 *> column j+1 of A, not counting the diagonal. Hence
197 *>
198 *> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
199 *> 1<=i<=j
200 *> and
201 *>
202 *> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
203 *> 1<=i< j
204 *>
205 *> Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTBSV if the
206 *> reciprocal of the largest M(j), j=1,..,n, is larger than
207 *> max(underflow, 1/overflow).
208 *>
209 *> The bound on x(j) is also used to determine when a step in the
210 *> columnwise method can be performed without fear of overflow. If
211 *> the computed bound is greater than a large constant, x is scaled to
212 *> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
213 *> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
214 *>
215 *> Similarly, a row-wise scheme is used to solve A**T*x = b. The basic
216 *> algorithm for A upper triangular is
217 *>
218 *> for j = 1, ..., n
219 *> x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
220 *> end
221 *>
222 *> We simultaneously compute two bounds
223 *> G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
224 *> M(j) = bound on x(i), 1<=i<=j
225 *>
226 *> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
227 *> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
228 *> Then the bound on x(j) is
229 *>
230 *> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
231 *>
232 *> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
233 *> 1<=i<=j
234 *>
235 *> and we can safely call DTBSV if 1/M(n) and 1/G(n) are both greater
236 *> than max(underflow, 1/overflow).
237 *> \endverbatim
238 *>
239 * =====================================================================
240  SUBROUTINE dlatbs( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
241  $ SCALE, CNORM, INFO )
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  DOUBLE PRECISION SCALE
251 * ..
252 * .. Array Arguments ..
253  DOUBLE PRECISION AB( LDAB, * ), CNORM( * ), X( * )
254 * ..
255 *
256 * =====================================================================
257 *
258 * .. Parameters ..
259  DOUBLE PRECISION ZERO, HALF, ONE
260  parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
261 * ..
262 * .. Local Scalars ..
263  LOGICAL NOTRAN, NOUNIT, UPPER
264  INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
265  DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
266  $ tmax, tscal, uscal, xbnd, xj, xmax
267 * ..
268 * .. External Functions ..
269  LOGICAL LSAME
270  INTEGER IDAMAX
271  DOUBLE PRECISION DASUM, DDOT, DLAMCH
272  EXTERNAL lsame, idamax, dasum, ddot, dlamch
273 * ..
274 * .. External Subroutines ..
275  EXTERNAL daxpy, dscal, dtbsv, 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( 'DLATBS', -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 = dlamch( 'Safe minimum' ) / dlamch( '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 ) = dasum( 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 ) = dasum( 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 = idamax( 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 dscal( 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 DTBSV can be used.
363 *
364  j = idamax( 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 dtbsv( 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 dscal( n, scale, x, 1 )
531  xmax = bignum
532  END IF
533 *
534  IF( notran ) THEN
535 *
536 * Solve A * x = b
537 *
538  DO 110 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 100
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 dscal( 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 dscal( 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  100 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 dscal( 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 dscal( 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 daxpy( jlen, -x( j )*tscal,
636  $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
637  i = idamax( 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 daxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
649  $ x( j+1 ), 1 )
650  i = j + idamax( n-j, x( j+1 ), 1 )
651  xmax = abs( x( i ) )
652  END IF
653  110 CONTINUE
654 *
655  ELSE
656 *
657 * Solve A**T * x = b
658 *
659  DO 160 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 dscal( 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 DDOT to perform the dot product.
697 *
698  IF( upper ) THEN
699  jlen = min( kd, j-1 )
700  sumj = ddot( 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 = ddot( 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 120 i = 1, jlen
714  sumj = sumj + ( ab( kd+i-jlen, j )*uscal )*
715  $ x( j-jlen-1+i )
716  120 CONTINUE
717  ELSE
718  jlen = min( kd, n-j )
719  DO 130 i = 1, jlen
720  sumj = sumj + ( ab( i+1, j )*uscal )*x( j+i )
721  130 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 150
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 dscal( 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 dscal( 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 140 i = 1, n
779  x( i ) = zero
780  140 CONTINUE
781  x( j ) = one
782  scale = zero
783  xmax = zero
784  END IF
785  150 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  160 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 dscal( n, one / tscal, cnorm, 1 )
803  END IF
804 *
805  RETURN
806 *
807 * End of DLATBS
808 *
809  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dtbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
DTBSV
Definition: dtbsv.f:189
subroutine dlatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
DLATBS solves a triangular banded system of equations.
Definition: dlatbs.f:242