LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
clatbs.f
Go to the documentation of this file.
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 *> \date September 2012
167 *
168 *> \ingroup complexOTHERauxiliary
169 *
170 *> \par Further Details:
171 * =====================
172 *>
173 *> \verbatim
174 *>
175 *> A rough bound on x is computed; if that is less than overflow, CTBSV
176 *> is called, otherwise, specific code is used which checks for possible
177 *> overflow or divide-by-zero at every operation.
178 *>
179 *> A columnwise scheme is used for solving A*x = b. The basic algorithm
180 *> if A is lower triangular is
181 *>
182 *> x[1:n] := b[1:n]
183 *> for j = 1, ..., n
184 *> x(j) := x(j) / A(j,j)
185 *> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
186 *> end
187 *>
188 *> Define bounds on the components of x after j iterations of the loop:
189 *> M(j) = bound on x[1:j]
190 *> G(j) = bound on x[j+1:n]
191 *> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
192 *>
193 *> Then for iteration j+1 we have
194 *> M(j+1) <= G(j) / | A(j+1,j+1) |
195 *> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
196 *> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
197 *>
198 *> where CNORM(j+1) is greater than or equal to the infinity-norm of
199 *> column j+1 of A, not counting the diagonal. Hence
200 *>
201 *> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
202 *> 1<=i<=j
203 *> and
204 *>
205 *> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
206 *> 1<=i< j
207 *>
208 *> Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTBSV if the
209 *> reciprocal of the largest M(j), j=1,..,n, is larger than
210 *> max(underflow, 1/overflow).
211 *>
212 *> The bound on x(j) is also used to determine when a step in the
213 *> columnwise method can be performed without fear of overflow. If
214 *> the computed bound is greater than a large constant, x is scaled to
215 *> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
216 *> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
217 *>
218 *> Similarly, a row-wise scheme is used to solve A**T *x = b or
219 *> A**H *x = b. The basic algorithm for A upper triangular is
220 *>
221 *> for j = 1, ..., n
222 *> x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
223 *> end
224 *>
225 *> We simultaneously compute two bounds
226 *> G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
227 *> M(j) = bound on x(i), 1<=i<=j
228 *>
229 *> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
230 *> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
231 *> Then the bound on x(j) is
232 *>
233 *> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
234 *>
235 *> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
236 *> 1<=i<=j
237 *>
238 *> and we can safely call CTBSV if 1/M(n) and 1/G(n) are both greater
239 *> than max(underflow, 1/overflow).
240 *> \endverbatim
241 *>
242 * =====================================================================
243  SUBROUTINE clatbs( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
244  $ scale, cnorm, info )
245 *
246 * -- LAPACK auxiliary routine (version 3.4.2) --
247 * -- LAPACK is a software package provided by Univ. of Tennessee, --
248 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
249 * September 2012
250 *
251 * .. Scalar Arguments ..
252  CHARACTER DIAG, NORMIN, TRANS, UPLO
253  INTEGER INFO, KD, LDAB, N
254  REAL SCALE
255 * ..
256 * .. Array Arguments ..
257  REAL CNORM( * )
258  COMPLEX AB( ldab, * ), X( * )
259 * ..
260 *
261 * =====================================================================
262 *
263 * .. Parameters ..
264  REAL ZERO, HALF, ONE, TWO
265  parameter ( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
266  $ two = 2.0e+0 )
267 * ..
268 * .. Local Scalars ..
269  LOGICAL NOTRAN, NOUNIT, UPPER
270  INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
271  REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
272  $ xbnd, xj, xmax
273  COMPLEX CSUMJ, TJJS, USCAL, ZDUM
274 * ..
275 * .. External Functions ..
276  LOGICAL LSAME
277  INTEGER ICAMAX, ISAMAX
278  REAL SCASUM, SLAMCH
279  COMPLEX CDOTC, CDOTU, CLADIV
280  EXTERNAL lsame, icamax, isamax, scasum, slamch, cdotc,
281  $ cdotu, cladiv
282 * ..
283 * .. External Subroutines ..
284  EXTERNAL caxpy, csscal, ctbsv, slabad, sscal, xerbla
285 * ..
286 * .. Intrinsic Functions ..
287  INTRINSIC abs, aimag, cmplx, conjg, max, min, real
288 * ..
289 * .. Statement Functions ..
290  REAL CABS1, CABS2
291 * ..
292 * .. Statement Function definitions ..
293  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( AIMAG( zdum ) )
294  cabs2( zdum ) = abs( REAL( ZDUM ) / 2. ) +
295  $ abs( aimag( zdum ) / 2. )
296 * ..
297 * .. Executable Statements ..
298 *
299  info = 0
300  upper = lsame( uplo, 'U' )
301  notran = lsame( trans, 'N' )
302  nounit = lsame( diag, 'N' )
303 *
304 * Test the input parameters.
305 *
306  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
307  info = -1
308  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
309  $ lsame( trans, 'C' ) ) THEN
310  info = -2
311  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
312  info = -3
313  ELSE IF( .NOT.lsame( normin, 'Y' ) .AND. .NOT.
314  $ lsame( normin, 'N' ) ) THEN
315  info = -4
316  ELSE IF( n.LT.0 ) THEN
317  info = -5
318  ELSE IF( kd.LT.0 ) THEN
319  info = -6
320  ELSE IF( ldab.LT.kd+1 ) THEN
321  info = -8
322  END IF
323  IF( info.NE.0 ) THEN
324  CALL xerbla( 'CLATBS', -info )
325  RETURN
326  END IF
327 *
328 * Quick return if possible
329 *
330  IF( n.EQ.0 )
331  $ RETURN
332 *
333 * Determine machine dependent parameters to control overflow.
334 *
335  smlnum = slamch( 'Safe minimum' )
336  bignum = one / smlnum
337  CALL slabad( smlnum, bignum )
338  smlnum = smlnum / slamch( 'Precision' )
339  bignum = one / smlnum
340  scale = one
341 *
342  IF( lsame( normin, 'N' ) ) THEN
343 *
344 * Compute the 1-norm of each column, not including the diagonal.
345 *
346  IF( upper ) THEN
347 *
348 * A is upper triangular.
349 *
350  DO 10 j = 1, n
351  jlen = min( kd, j-1 )
352  cnorm( j ) = scasum( jlen, ab( kd+1-jlen, j ), 1 )
353  10 CONTINUE
354  ELSE
355 *
356 * A is lower triangular.
357 *
358  DO 20 j = 1, n
359  jlen = min( kd, n-j )
360  IF( jlen.GT.0 ) THEN
361  cnorm( j ) = scasum( jlen, ab( 2, j ), 1 )
362  ELSE
363  cnorm( j ) = zero
364  END IF
365  20 CONTINUE
366  END IF
367  END IF
368 *
369 * Scale the column norms by TSCAL if the maximum element in CNORM is
370 * greater than BIGNUM/2.
371 *
372  imax = isamax( n, cnorm, 1 )
373  tmax = cnorm( imax )
374  IF( tmax.LE.bignum*half ) THEN
375  tscal = one
376  ELSE
377  tscal = half / ( smlnum*tmax )
378  CALL sscal( n, tscal, cnorm, 1 )
379  END IF
380 *
381 * Compute a bound on the computed solution vector to see if the
382 * Level 2 BLAS routine CTBSV can be used.
383 *
384  xmax = zero
385  DO 30 j = 1, n
386  xmax = max( xmax, cabs2( x( j ) ) )
387  30 CONTINUE
388  xbnd = xmax
389  IF( notran ) THEN
390 *
391 * Compute the growth in A * x = b.
392 *
393  IF( upper ) THEN
394  jfirst = n
395  jlast = 1
396  jinc = -1
397  maind = kd + 1
398  ELSE
399  jfirst = 1
400  jlast = n
401  jinc = 1
402  maind = 1
403  END IF
404 *
405  IF( tscal.NE.one ) THEN
406  grow = zero
407  GO TO 60
408  END IF
409 *
410  IF( nounit ) THEN
411 *
412 * A is non-unit triangular.
413 *
414 * Compute GROW = 1/G(j) and XBND = 1/M(j).
415 * Initially, G(0) = max{x(i), i=1,...,n}.
416 *
417  grow = half / max( xbnd, smlnum )
418  xbnd = grow
419  DO 40 j = jfirst, jlast, jinc
420 *
421 * Exit the loop if the growth factor is too small.
422 *
423  IF( grow.LE.smlnum )
424  $ GO TO 60
425 *
426  tjjs = ab( maind, j )
427  tjj = cabs1( tjjs )
428 *
429  IF( tjj.GE.smlnum ) THEN
430 *
431 * M(j) = G(j-1) / abs(A(j,j))
432 *
433  xbnd = min( xbnd, min( one, tjj )*grow )
434  ELSE
435 *
436 * M(j) could overflow, set XBND to 0.
437 *
438  xbnd = zero
439  END IF
440 *
441  IF( tjj+cnorm( j ).GE.smlnum ) THEN
442 *
443 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
444 *
445  grow = grow*( tjj / ( tjj+cnorm( j ) ) )
446  ELSE
447 *
448 * G(j) could overflow, set GROW to 0.
449 *
450  grow = zero
451  END IF
452  40 CONTINUE
453  grow = xbnd
454  ELSE
455 *
456 * A is unit triangular.
457 *
458 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
459 *
460  grow = min( one, half / max( xbnd, smlnum ) )
461  DO 50 j = jfirst, jlast, jinc
462 *
463 * Exit the loop if the growth factor is too small.
464 *
465  IF( grow.LE.smlnum )
466  $ GO TO 60
467 *
468 * G(j) = G(j-1)*( 1 + CNORM(j) )
469 *
470  grow = grow*( one / ( one+cnorm( j ) ) )
471  50 CONTINUE
472  END IF
473  60 CONTINUE
474 *
475  ELSE
476 *
477 * Compute the growth in A**T * x = b or A**H * x = b.
478 *
479  IF( upper ) THEN
480  jfirst = 1
481  jlast = n
482  jinc = 1
483  maind = kd + 1
484  ELSE
485  jfirst = n
486  jlast = 1
487  jinc = -1
488  maind = 1
489  END IF
490 *
491  IF( tscal.NE.one ) THEN
492  grow = zero
493  GO TO 90
494  END IF
495 *
496  IF( nounit ) THEN
497 *
498 * A is non-unit triangular.
499 *
500 * Compute GROW = 1/G(j) and XBND = 1/M(j).
501 * Initially, M(0) = max{x(i), i=1,...,n}.
502 *
503  grow = half / max( xbnd, smlnum )
504  xbnd = grow
505  DO 70 j = jfirst, jlast, jinc
506 *
507 * Exit the loop if the growth factor is too small.
508 *
509  IF( grow.LE.smlnum )
510  $ GO TO 90
511 *
512 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
513 *
514  xj = one + cnorm( j )
515  grow = min( grow, xbnd / xj )
516 *
517  tjjs = ab( maind, j )
518  tjj = cabs1( tjjs )
519 *
520  IF( tjj.GE.smlnum ) THEN
521 *
522 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
523 *
524  IF( xj.GT.tjj )
525  $ xbnd = xbnd*( tjj / xj )
526  ELSE
527 *
528 * M(j) could overflow, set XBND to 0.
529 *
530  xbnd = zero
531  END IF
532  70 CONTINUE
533  grow = min( grow, xbnd )
534  ELSE
535 *
536 * A is unit triangular.
537 *
538 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
539 *
540  grow = min( one, half / max( xbnd, smlnum ) )
541  DO 80 j = jfirst, jlast, jinc
542 *
543 * Exit the loop if the growth factor is too small.
544 *
545  IF( grow.LE.smlnum )
546  $ GO TO 90
547 *
548 * G(j) = ( 1 + CNORM(j) )*G(j-1)
549 *
550  xj = one + cnorm( j )
551  grow = grow / xj
552  80 CONTINUE
553  END IF
554  90 CONTINUE
555  END IF
556 *
557  IF( ( grow*tscal ).GT.smlnum ) THEN
558 *
559 * Use the Level 2 BLAS solve if the reciprocal of the bound on
560 * elements of X is not too small.
561 *
562  CALL ctbsv( uplo, trans, diag, n, kd, ab, ldab, x, 1 )
563  ELSE
564 *
565 * Use a Level 1 BLAS solve, scaling intermediate results.
566 *
567  IF( xmax.GT.bignum*half ) THEN
568 *
569 * Scale X so that its components are less than or equal to
570 * BIGNUM in absolute value.
571 *
572  scale = ( bignum*half ) / xmax
573  CALL csscal( n, scale, x, 1 )
574  xmax = bignum
575  ELSE
576  xmax = xmax*two
577  END IF
578 *
579  IF( notran ) THEN
580 *
581 * Solve A * x = b
582 *
583  DO 110 j = jfirst, jlast, jinc
584 *
585 * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
586 *
587  xj = cabs1( x( j ) )
588  IF( nounit ) THEN
589  tjjs = ab( maind, j )*tscal
590  ELSE
591  tjjs = tscal
592  IF( tscal.EQ.one )
593  $ GO TO 105
594  END IF
595  tjj = cabs1( tjjs )
596  IF( tjj.GT.smlnum ) THEN
597 *
598 * abs(A(j,j)) > SMLNUM:
599 *
600  IF( tjj.LT.one ) THEN
601  IF( xj.GT.tjj*bignum ) THEN
602 *
603 * Scale x by 1/b(j).
604 *
605  rec = one / xj
606  CALL csscal( n, rec, x, 1 )
607  scale = scale*rec
608  xmax = xmax*rec
609  END IF
610  END IF
611  x( j ) = cladiv( x( j ), tjjs )
612  xj = cabs1( x( j ) )
613  ELSE IF( tjj.GT.zero ) THEN
614 *
615 * 0 < abs(A(j,j)) <= SMLNUM:
616 *
617  IF( xj.GT.tjj*bignum ) THEN
618 *
619 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
620 * to avoid overflow when dividing by A(j,j).
621 *
622  rec = ( tjj*bignum ) / xj
623  IF( cnorm( j ).GT.one ) THEN
624 *
625 * Scale by 1/CNORM(j) to avoid overflow when
626 * multiplying x(j) times column j.
627 *
628  rec = rec / cnorm( j )
629  END IF
630  CALL csscal( n, rec, x, 1 )
631  scale = scale*rec
632  xmax = xmax*rec
633  END IF
634  x( j ) = cladiv( x( j ), tjjs )
635  xj = cabs1( x( j ) )
636  ELSE
637 *
638 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
639 * scale = 0, and compute a solution to A*x = 0.
640 *
641  DO 100 i = 1, n
642  x( i ) = zero
643  100 CONTINUE
644  x( j ) = one
645  xj = one
646  scale = zero
647  xmax = zero
648  END IF
649  105 CONTINUE
650 *
651 * Scale x if necessary to avoid overflow when adding a
652 * multiple of column j of A.
653 *
654  IF( xj.GT.one ) THEN
655  rec = one / xj
656  IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
657 *
658 * Scale x by 1/(2*abs(x(j))).
659 *
660  rec = rec*half
661  CALL csscal( n, rec, x, 1 )
662  scale = scale*rec
663  END IF
664  ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
665 *
666 * Scale x by 1/2.
667 *
668  CALL csscal( n, half, x, 1 )
669  scale = scale*half
670  END IF
671 *
672  IF( upper ) THEN
673  IF( j.GT.1 ) THEN
674 *
675 * Compute the update
676 * x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
677 * x(j)* A(max(1,j-kd):j-1,j)
678 *
679  jlen = min( kd, j-1 )
680  CALL caxpy( jlen, -x( j )*tscal,
681  $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
682  i = icamax( j-1, x, 1 )
683  xmax = cabs1( x( i ) )
684  END IF
685  ELSE IF( j.LT.n ) THEN
686 *
687 * Compute the update
688 * x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
689 * x(j) * A(j+1:min(j+kd,n),j)
690 *
691  jlen = min( kd, n-j )
692  IF( jlen.GT.0 )
693  $ CALL caxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
694  $ x( j+1 ), 1 )
695  i = j + icamax( n-j, x( j+1 ), 1 )
696  xmax = cabs1( x( i ) )
697  END IF
698  110 CONTINUE
699 *
700  ELSE IF( lsame( trans, 'T' ) ) THEN
701 *
702 * Solve A**T * x = b
703 *
704  DO 150 j = jfirst, jlast, jinc
705 *
706 * Compute x(j) = b(j) - sum A(k,j)*x(k).
707 * k<>j
708 *
709  xj = cabs1( x( j ) )
710  uscal = tscal
711  rec = one / max( xmax, one )
712  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
713 *
714 * If x(j) could overflow, scale x by 1/(2*XMAX).
715 *
716  rec = rec*half
717  IF( nounit ) THEN
718  tjjs = ab( maind, j )*tscal
719  ELSE
720  tjjs = tscal
721  END IF
722  tjj = cabs1( tjjs )
723  IF( tjj.GT.one ) THEN
724 *
725 * Divide by A(j,j) when scaling x if A(j,j) > 1.
726 *
727  rec = min( one, rec*tjj )
728  uscal = cladiv( uscal, tjjs )
729  END IF
730  IF( rec.LT.one ) THEN
731  CALL csscal( n, rec, x, 1 )
732  scale = scale*rec
733  xmax = xmax*rec
734  END IF
735  END IF
736 *
737  csumj = zero
738  IF( uscal.EQ.cmplx( one ) ) THEN
739 *
740 * If the scaling needed for A in the dot product is 1,
741 * call CDOTU to perform the dot product.
742 *
743  IF( upper ) THEN
744  jlen = min( kd, j-1 )
745  csumj = cdotu( jlen, ab( kd+1-jlen, j ), 1,
746  $ x( j-jlen ), 1 )
747  ELSE
748  jlen = min( kd, n-j )
749  IF( jlen.GT.1 )
750  $ csumj = cdotu( jlen, ab( 2, j ), 1, x( j+1 ),
751  $ 1 )
752  END IF
753  ELSE
754 *
755 * Otherwise, use in-line code for the dot product.
756 *
757  IF( upper ) THEN
758  jlen = min( kd, j-1 )
759  DO 120 i = 1, jlen
760  csumj = csumj + ( ab( kd+i-jlen, j )*uscal )*
761  $ x( j-jlen-1+i )
762  120 CONTINUE
763  ELSE
764  jlen = min( kd, n-j )
765  DO 130 i = 1, jlen
766  csumj = csumj + ( ab( i+1, j )*uscal )*x( j+i )
767  130 CONTINUE
768  END IF
769  END IF
770 *
771  IF( uscal.EQ.cmplx( tscal ) ) THEN
772 *
773 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
774 * was not used to scale the dotproduct.
775 *
776  x( j ) = x( j ) - csumj
777  xj = cabs1( x( j ) )
778  IF( nounit ) THEN
779 *
780 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
781 *
782  tjjs = ab( maind, j )*tscal
783  ELSE
784  tjjs = tscal
785  IF( tscal.EQ.one )
786  $ GO TO 145
787  END IF
788  tjj = cabs1( tjjs )
789  IF( tjj.GT.smlnum ) THEN
790 *
791 * abs(A(j,j)) > SMLNUM:
792 *
793  IF( tjj.LT.one ) THEN
794  IF( xj.GT.tjj*bignum ) THEN
795 *
796 * Scale X by 1/abs(x(j)).
797 *
798  rec = one / xj
799  CALL csscal( n, rec, x, 1 )
800  scale = scale*rec
801  xmax = xmax*rec
802  END IF
803  END IF
804  x( j ) = cladiv( x( j ), tjjs )
805  ELSE IF( tjj.GT.zero ) THEN
806 *
807 * 0 < abs(A(j,j)) <= SMLNUM:
808 *
809  IF( xj.GT.tjj*bignum ) THEN
810 *
811 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
812 *
813  rec = ( tjj*bignum ) / xj
814  CALL csscal( n, rec, x, 1 )
815  scale = scale*rec
816  xmax = xmax*rec
817  END IF
818  x( j ) = cladiv( x( j ), tjjs )
819  ELSE
820 *
821 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
822 * scale = 0 and compute a solution to A**T *x = 0.
823 *
824  DO 140 i = 1, n
825  x( i ) = zero
826  140 CONTINUE
827  x( j ) = one
828  scale = zero
829  xmax = zero
830  END IF
831  145 CONTINUE
832  ELSE
833 *
834 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
835 * product has already been divided by 1/A(j,j).
836 *
837  x( j ) = cladiv( x( j ), tjjs ) - csumj
838  END IF
839  xmax = max( xmax, cabs1( x( j ) ) )
840  150 CONTINUE
841 *
842  ELSE
843 *
844 * Solve A**H * x = b
845 *
846  DO 190 j = jfirst, jlast, jinc
847 *
848 * Compute x(j) = b(j) - sum A(k,j)*x(k).
849 * k<>j
850 *
851  xj = cabs1( x( j ) )
852  uscal = tscal
853  rec = one / max( xmax, one )
854  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
855 *
856 * If x(j) could overflow, scale x by 1/(2*XMAX).
857 *
858  rec = rec*half
859  IF( nounit ) THEN
860  tjjs = conjg( ab( maind, j ) )*tscal
861  ELSE
862  tjjs = tscal
863  END IF
864  tjj = cabs1( tjjs )
865  IF( tjj.GT.one ) THEN
866 *
867 * Divide by A(j,j) when scaling x if A(j,j) > 1.
868 *
869  rec = min( one, rec*tjj )
870  uscal = cladiv( uscal, tjjs )
871  END IF
872  IF( rec.LT.one ) THEN
873  CALL csscal( n, rec, x, 1 )
874  scale = scale*rec
875  xmax = xmax*rec
876  END IF
877  END IF
878 *
879  csumj = zero
880  IF( uscal.EQ.cmplx( one ) ) THEN
881 *
882 * If the scaling needed for A in the dot product is 1,
883 * call CDOTC to perform the dot product.
884 *
885  IF( upper ) THEN
886  jlen = min( kd, j-1 )
887  csumj = cdotc( jlen, ab( kd+1-jlen, j ), 1,
888  $ x( j-jlen ), 1 )
889  ELSE
890  jlen = min( kd, n-j )
891  IF( jlen.GT.1 )
892  $ csumj = cdotc( jlen, ab( 2, j ), 1, x( j+1 ),
893  $ 1 )
894  END IF
895  ELSE
896 *
897 * Otherwise, use in-line code for the dot product.
898 *
899  IF( upper ) THEN
900  jlen = min( kd, j-1 )
901  DO 160 i = 1, jlen
902  csumj = csumj + ( conjg( ab( kd+i-jlen, j ) )*
903  $ uscal )*x( j-jlen-1+i )
904  160 CONTINUE
905  ELSE
906  jlen = min( kd, n-j )
907  DO 170 i = 1, jlen
908  csumj = csumj + ( conjg( ab( i+1, j ) )*uscal )*
909  $ x( j+i )
910  170 CONTINUE
911  END IF
912  END IF
913 *
914  IF( uscal.EQ.cmplx( tscal ) ) THEN
915 *
916 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
917 * was not used to scale the dotproduct.
918 *
919  x( j ) = x( j ) - csumj
920  xj = cabs1( x( j ) )
921  IF( nounit ) THEN
922 *
923 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
924 *
925  tjjs = conjg( ab( maind, j ) )*tscal
926  ELSE
927  tjjs = tscal
928  IF( tscal.EQ.one )
929  $ GO TO 185
930  END IF
931  tjj = cabs1( tjjs )
932  IF( tjj.GT.smlnum ) THEN
933 *
934 * abs(A(j,j)) > SMLNUM:
935 *
936  IF( tjj.LT.one ) THEN
937  IF( xj.GT.tjj*bignum ) THEN
938 *
939 * Scale X by 1/abs(x(j)).
940 *
941  rec = one / xj
942  CALL csscal( n, rec, x, 1 )
943  scale = scale*rec
944  xmax = xmax*rec
945  END IF
946  END IF
947  x( j ) = cladiv( x( j ), tjjs )
948  ELSE IF( tjj.GT.zero ) THEN
949 *
950 * 0 < abs(A(j,j)) <= SMLNUM:
951 *
952  IF( xj.GT.tjj*bignum ) THEN
953 *
954 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
955 *
956  rec = ( tjj*bignum ) / xj
957  CALL csscal( n, rec, x, 1 )
958  scale = scale*rec
959  xmax = xmax*rec
960  END IF
961  x( j ) = cladiv( x( j ), tjjs )
962  ELSE
963 *
964 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
965 * scale = 0 and compute a solution to A**H *x = 0.
966 *
967  DO 180 i = 1, n
968  x( i ) = zero
969  180 CONTINUE
970  x( j ) = one
971  scale = zero
972  xmax = zero
973  END IF
974  185 CONTINUE
975  ELSE
976 *
977 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
978 * product has already been divided by 1/A(j,j).
979 *
980  x( j ) = cladiv( x( j ), tjjs ) - csumj
981  END IF
982  xmax = max( xmax, cabs1( x( j ) ) )
983  190 CONTINUE
984  END IF
985  scale = scale / tscal
986  END IF
987 *
988 * Scale the column norms by 1/TSCAL for return.
989 *
990  IF( tscal.NE.one ) THEN
991  CALL sscal( n, one / tscal, cnorm, 1 )
992  END IF
993 *
994  RETURN
995 *
996 * End of CLATBS
997 *
998  END
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
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:245
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
ctbsv
subroutine ctbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
CTBSV
Definition: ctbsv.f:191
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54