LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
clatbs.f
1 *> \brief \b CLATBS 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 CLATBS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatbs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatbs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatbs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLATBS( 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 CNORM( * )
31 * COMPLEX AB( LDAB, * ), X( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CLATBS solves one of the triangular systems
41 *>
42 *> A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
43 *>
44 *> with scaling to prevent overflow, where A is an upper or lower
45 *> triangular band matrix. Here A**T denotes the transpose of A, x and b
46 *> are n-element vectors, and s is a scaling factor, usually less than
47 *> or equal to 1, chosen so that the components of x will be less than
48 *> the overflow threshold. If the unscaled problem will not cause
49 *> overflow, the Level 2 BLAS routine CTBSV is called. If the matrix A
50 *> is singular (A(j,j) = 0 for some j), then s is set to 0 and a
51 *> non-trivial solution to A*x = 0 is returned.
52 *> \endverbatim
53 *
54 * Arguments:
55 * ==========
56 *
57 *> \param[in] UPLO
58 *> \verbatim
59 *> UPLO is CHARACTER*1
60 *> Specifies whether the matrix A is upper or lower triangular.
61 *> = 'U': Upper triangular
62 *> = 'L': Lower triangular
63 *> \endverbatim
64 *>
65 *> \param[in] TRANS
66 *> \verbatim
67 *> TRANS is CHARACTER*1
68 *> Specifies the operation applied to A.
69 *> = 'N': Solve A * x = s*b (No transpose)
70 *> = 'T': Solve A**T * x = s*b (Transpose)
71 *> = 'C': Solve A**H * x = s*b (Conjugate transpose)
72 *> \endverbatim
73 *>
74 *> \param[in] DIAG
75 *> \verbatim
76 *> DIAG is CHARACTER*1
77 *> Specifies whether or not the matrix A is unit triangular.
78 *> = 'N': Non-unit triangular
79 *> = 'U': Unit triangular
80 *> \endverbatim
81 *>
82 *> \param[in] NORMIN
83 *> \verbatim
84 *> NORMIN is CHARACTER*1
85 *> Specifies whether CNORM has been set or not.
86 *> = 'Y': CNORM contains the column norms on entry
87 *> = 'N': CNORM is not set on entry. On exit, the norms will
88 *> be computed and stored in CNORM.
89 *> \endverbatim
90 *>
91 *> \param[in] N
92 *> \verbatim
93 *> N is INTEGER
94 *> The order of the matrix A. N >= 0.
95 *> \endverbatim
96 *>
97 *> \param[in] KD
98 *> \verbatim
99 *> KD is INTEGER
100 *> The number of subdiagonals or superdiagonals in the
101 *> triangular matrix A. KD >= 0.
102 *> \endverbatim
103 *>
104 *> \param[in] AB
105 *> \verbatim
106 *> AB is COMPLEX array, dimension (LDAB,N)
107 *> The upper or lower triangular band matrix A, stored in the
108 *> first KD+1 rows of the array. The j-th column of A is stored
109 *> in the j-th column of the array AB as follows:
110 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
111 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
112 *> \endverbatim
113 *>
114 *> \param[in] LDAB
115 *> \verbatim
116 *> LDAB is INTEGER
117 *> The leading dimension of the array AB. LDAB >= KD+1.
118 *> \endverbatim
119 *>
120 *> \param[in,out] X
121 *> \verbatim
122 *> X is COMPLEX array, dimension (N)
123 *> On entry, the right hand side b of the triangular system.
124 *> On exit, X is overwritten by the solution vector x.
125 *> \endverbatim
126 *>
127 *> \param[out] SCALE
128 *> \verbatim
129 *> SCALE is REAL
130 *> The scaling factor s for the triangular system
131 *> A * x = s*b, A**T * x = s*b, or A**H * x = s*b.
132 *> If SCALE = 0, the matrix A is singular or badly scaled, and
133 *> the vector x is an exact or approximate solution to A*x = 0.
134 *> \endverbatim
135 *>
136 *> \param[in,out] CNORM
137 *> \verbatim
138 *> CNORM is REAL array, dimension (N)
139 *>
140 *> If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
141 *> contains the norm of the off-diagonal part of the j-th column
142 *> of A. If TRANS = 'N', CNORM(j) must be greater than or equal
143 *> to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
144 *> must be greater than or equal to the 1-norm.
145 *>
146 *> If NORMIN = 'N', CNORM is an output argument and CNORM(j)
147 *> returns the 1-norm of the offdiagonal part of the j-th column
148 *> of A.
149 *> \endverbatim
150 *>
151 *> \param[out] INFO
152 *> \verbatim
153 *> INFO is INTEGER
154 *> = 0: successful exit
155 *> < 0: if INFO = -k, the k-th argument had an illegal value
156 *> \endverbatim
157 *
158 * Authors:
159 * ========
160 *
161 *> \author Univ. of Tennessee
162 *> \author Univ. of California Berkeley
163 *> \author Univ. of Colorado Denver
164 *> \author NAG Ltd.
165 *
166 *> \ingroup latbs
167 *
168 *> \par Further Details:
169 * =====================
170 *>
171 *> \verbatim
172 *>
173 *> A rough bound on x is computed; if that is less than overflow, CTBSV
174 *> is called, otherwise, specific code is used which checks for possible
175 *> overflow or divide-by-zero at every operation.
176 *>
177 *> A columnwise scheme is used for solving A*x = b. The basic algorithm
178 *> if A is lower triangular is
179 *>
180 *> x[1:n] := b[1:n]
181 *> for j = 1, ..., n
182 *> x(j) := x(j) / A(j,j)
183 *> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
184 *> end
185 *>
186 *> Define bounds on the components of x after j iterations of the loop:
187 *> M(j) = bound on x[1:j]
188 *> G(j) = bound on x[j+1:n]
189 *> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
190 *>
191 *> Then for iteration j+1 we have
192 *> M(j+1) <= G(j) / | A(j+1,j+1) |
193 *> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
194 *> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
195 *>
196 *> where CNORM(j+1) is greater than or equal to the infinity-norm of
197 *> column j+1 of A, not counting the diagonal. Hence
198 *>
199 *> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
200 *> 1<=i<=j
201 *> and
202 *>
203 *> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
204 *> 1<=i< j
205 *>
206 *> Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTBSV if the
207 *> reciprocal of the largest M(j), j=1,..,n, is larger than
208 *> max(underflow, 1/overflow).
209 *>
210 *> The bound on x(j) is also used to determine when a step in the
211 *> columnwise method can be performed without fear of overflow. If
212 *> the computed bound is greater than a large constant, x is scaled to
213 *> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
214 *> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
215 *>
216 *> Similarly, a row-wise scheme is used to solve A**T *x = b or
217 *> A**H *x = b. The basic algorithm for A upper triangular is
218 *>
219 *> for j = 1, ..., n
220 *> x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
221 *> end
222 *>
223 *> We simultaneously compute two bounds
224 *> G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
225 *> M(j) = bound on x(i), 1<=i<=j
226 *>
227 *> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
228 *> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
229 *> Then the bound on x(j) is
230 *>
231 *> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
232 *>
233 *> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
234 *> 1<=i<=j
235 *>
236 *> and we can safely call CTBSV if 1/M(n) and 1/G(n) are both greater
237 *> than max(underflow, 1/overflow).
238 *> \endverbatim
239 *>
240 * =====================================================================
241  SUBROUTINE clatbs( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
242  $ SCALE, CNORM, INFO )
243 *
244 * -- LAPACK auxiliary routine --
245 * -- LAPACK is a software package provided by Univ. of Tennessee, --
246 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
247 *
248 * .. Scalar Arguments ..
249  CHARACTER DIAG, NORMIN, TRANS, UPLO
250  INTEGER INFO, KD, LDAB, N
251  REAL SCALE
252 * ..
253 * .. Array Arguments ..
254  REAL CNORM( * )
255  COMPLEX AB( ldab, * ), X( * )
256 * ..
257 *
258 * =====================================================================
259 *
260 * .. Parameters ..
261  REAL ZERO, HALF, ONE, TWO
262  parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
263  $ two = 2.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, TJJ, TMAX, TSCAL,
269  $ xbnd, xj, xmax
270  COMPLEX CSUMJ, TJJS, USCAL, ZDUM
271 * ..
272 * .. External Functions ..
273  LOGICAL LSAME
274  INTEGER ICAMAX, ISAMAX
275  REAL SCASUM, SLAMCH
276  COMPLEX CDOTC, CDOTU, CLADIV
277  EXTERNAL lsame, icamax, isamax, scasum, slamch, 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 *
992  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
clatbs
subroutine clatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
CLATBS solves a triangular banded system of equations.
Definition: clatbs.f:243
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
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88