LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zlatbs.f
Go to the documentation of this file.
1 *> \brief \b ZLATBS 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 ZLATBS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlatbs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlatbs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlatbs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLATBS( 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 CNORM( * )
31 * COMPLEX*16 AB( LDAB, * ), X( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> ZLATBS 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 ZTBSV 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*16 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*16 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 complex16OTHERauxiliary
167 *
168 *> \par Further Details:
169 * =====================
170 *>
171 *> \verbatim
172 *>
173 *> A rough bound on x is computed; if that is less than overflow, ZTBSV
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 ZTBSV 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 ZTBSV if 1/M(n) and 1/G(n) are both greater
237 *> than max(underflow, 1/overflow).
238 *> \endverbatim
239 *>
240 * =====================================================================
241  SUBROUTINE zlatbs( 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  DOUBLE PRECISION SCALE
252 * ..
253 * .. Array Arguments ..
254  DOUBLE PRECISION CNORM( * )
255  COMPLEX*16 AB( LDAB, * ), X( * )
256 * ..
257 *
258 * =====================================================================
259 *
260 * .. Parameters ..
261  DOUBLE PRECISION ZERO, HALF, ONE, TWO
262  parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0,
263  $ two = 2.0d+0 )
264 * ..
265 * .. Local Scalars ..
266  LOGICAL NOTRAN, NOUNIT, UPPER
267  INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
268  DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
269  $ xbnd, xj, xmax
270  COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
271 * ..
272 * .. External Functions ..
273  LOGICAL LSAME
274  INTEGER IDAMAX, IZAMAX
275  DOUBLE PRECISION DLAMCH, DZASUM
276  COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
277  EXTERNAL lsame, idamax, izamax, dlamch, dzasum, zdotc,
278  $ zdotu, zladiv
279 * ..
280 * .. External Subroutines ..
281  EXTERNAL dscal, xerbla, zaxpy, zdscal, ztbsv, dlabad
282 * ..
283 * .. Intrinsic Functions ..
284  INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
285 * ..
286 * .. Statement Functions ..
287  DOUBLE PRECISION CABS1, CABS2
288 * ..
289 * .. Statement Function definitions ..
290  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
291  cabs2( zdum ) = abs( dble( zdum ) / 2.d0 ) +
292  $ abs( dimag( zdum ) / 2.d0 )
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( 'ZLATBS', -info )
322  RETURN
323  END IF
324 *
325 * Quick return if possible
326 *
327  IF( n.EQ.0 )
328  $ RETURN
329 *
330 * Determine machine dependent parameters to control overflow.
331 *
332  smlnum = dlamch( 'Safe minimum' )
333  bignum = one / smlnum
334  CALL dlabad( smlnum, bignum )
335  smlnum = smlnum / dlamch( 'Precision' )
336  bignum = one / smlnum
337  scale = one
338 *
339  IF( lsame( normin, 'N' ) ) THEN
340 *
341 * Compute the 1-norm of each column, not including the diagonal.
342 *
343  IF( upper ) THEN
344 *
345 * A is upper triangular.
346 *
347  DO 10 j = 1, n
348  jlen = min( kd, j-1 )
349  cnorm( j ) = dzasum( jlen, ab( kd+1-jlen, j ), 1 )
350  10 CONTINUE
351  ELSE
352 *
353 * A is lower triangular.
354 *
355  DO 20 j = 1, n
356  jlen = min( kd, n-j )
357  IF( jlen.GT.0 ) THEN
358  cnorm( j ) = dzasum( jlen, ab( 2, j ), 1 )
359  ELSE
360  cnorm( j ) = zero
361  END IF
362  20 CONTINUE
363  END IF
364  END IF
365 *
366 * Scale the column norms by TSCAL if the maximum element in CNORM is
367 * greater than BIGNUM/2.
368 *
369  imax = idamax( n, cnorm, 1 )
370  tmax = cnorm( imax )
371  IF( tmax.LE.bignum*half ) THEN
372  tscal = one
373  ELSE
374  tscal = half / ( smlnum*tmax )
375  CALL dscal( n, tscal, cnorm, 1 )
376  END IF
377 *
378 * Compute a bound on the computed solution vector to see if the
379 * Level 2 BLAS routine ZTBSV can be used.
380 *
381  xmax = zero
382  DO 30 j = 1, n
383  xmax = max( xmax, cabs2( x( j ) ) )
384  30 CONTINUE
385  xbnd = xmax
386  IF( notran ) THEN
387 *
388 * Compute the growth in A * x = b.
389 *
390  IF( upper ) THEN
391  jfirst = n
392  jlast = 1
393  jinc = -1
394  maind = kd + 1
395  ELSE
396  jfirst = 1
397  jlast = n
398  jinc = 1
399  maind = 1
400  END IF
401 *
402  IF( tscal.NE.one ) THEN
403  grow = zero
404  GO TO 60
405  END IF
406 *
407  IF( nounit ) THEN
408 *
409 * A is non-unit triangular.
410 *
411 * Compute GROW = 1/G(j) and XBND = 1/M(j).
412 * Initially, G(0) = max{x(i), i=1,...,n}.
413 *
414  grow = half / max( xbnd, smlnum )
415  xbnd = grow
416  DO 40 j = jfirst, jlast, jinc
417 *
418 * Exit the loop if the growth factor is too small.
419 *
420  IF( grow.LE.smlnum )
421  $ GO TO 60
422 *
423  tjjs = ab( maind, j )
424  tjj = cabs1( tjjs )
425 *
426  IF( tjj.GE.smlnum ) THEN
427 *
428 * M(j) = G(j-1) / abs(A(j,j))
429 *
430  xbnd = min( xbnd, min( one, tjj )*grow )
431  ELSE
432 *
433 * M(j) could overflow, set XBND to 0.
434 *
435  xbnd = zero
436  END IF
437 *
438  IF( tjj+cnorm( j ).GE.smlnum ) THEN
439 *
440 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
441 *
442  grow = grow*( tjj / ( tjj+cnorm( j ) ) )
443  ELSE
444 *
445 * G(j) could overflow, set GROW to 0.
446 *
447  grow = zero
448  END IF
449  40 CONTINUE
450  grow = xbnd
451  ELSE
452 *
453 * A is unit triangular.
454 *
455 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
456 *
457  grow = min( one, half / max( xbnd, smlnum ) )
458  DO 50 j = jfirst, jlast, jinc
459 *
460 * Exit the loop if the growth factor is too small.
461 *
462  IF( grow.LE.smlnum )
463  $ GO TO 60
464 *
465 * G(j) = G(j-1)*( 1 + CNORM(j) )
466 *
467  grow = grow*( one / ( one+cnorm( j ) ) )
468  50 CONTINUE
469  END IF
470  60 CONTINUE
471 *
472  ELSE
473 *
474 * Compute the growth in A**T * x = b or A**H * x = b.
475 *
476  IF( upper ) THEN
477  jfirst = 1
478  jlast = n
479  jinc = 1
480  maind = kd + 1
481  ELSE
482  jfirst = n
483  jlast = 1
484  jinc = -1
485  maind = 1
486  END IF
487 *
488  IF( tscal.NE.one ) THEN
489  grow = zero
490  GO TO 90
491  END IF
492 *
493  IF( nounit ) THEN
494 *
495 * A is non-unit triangular.
496 *
497 * Compute GROW = 1/G(j) and XBND = 1/M(j).
498 * Initially, M(0) = max{x(i), i=1,...,n}.
499 *
500  grow = half / max( xbnd, smlnum )
501  xbnd = grow
502  DO 70 j = jfirst, jlast, jinc
503 *
504 * Exit the loop if the growth factor is too small.
505 *
506  IF( grow.LE.smlnum )
507  $ GO TO 90
508 *
509 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
510 *
511  xj = one + cnorm( j )
512  grow = min( grow, xbnd / xj )
513 *
514  tjjs = ab( maind, j )
515  tjj = cabs1( tjjs )
516 *
517  IF( tjj.GE.smlnum ) THEN
518 *
519 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
520 *
521  IF( xj.GT.tjj )
522  $ xbnd = xbnd*( tjj / xj )
523  ELSE
524 *
525 * M(j) could overflow, set XBND to 0.
526 *
527  xbnd = zero
528  END IF
529  70 CONTINUE
530  grow = min( grow, xbnd )
531  ELSE
532 *
533 * A is unit triangular.
534 *
535 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
536 *
537  grow = min( one, half / max( xbnd, smlnum ) )
538  DO 80 j = jfirst, jlast, jinc
539 *
540 * Exit the loop if the growth factor is too small.
541 *
542  IF( grow.LE.smlnum )
543  $ GO TO 90
544 *
545 * G(j) = ( 1 + CNORM(j) )*G(j-1)
546 *
547  xj = one + cnorm( j )
548  grow = grow / xj
549  80 CONTINUE
550  END IF
551  90 CONTINUE
552  END IF
553 *
554  IF( ( grow*tscal ).GT.smlnum ) THEN
555 *
556 * Use the Level 2 BLAS solve if the reciprocal of the bound on
557 * elements of X is not too small.
558 *
559  CALL ztbsv( uplo, trans, diag, n, kd, ab, ldab, x, 1 )
560  ELSE
561 *
562 * Use a Level 1 BLAS solve, scaling intermediate results.
563 *
564  IF( xmax.GT.bignum*half ) THEN
565 *
566 * Scale X so that its components are less than or equal to
567 * BIGNUM in absolute value.
568 *
569  scale = ( bignum*half ) / xmax
570  CALL zdscal( n, scale, x, 1 )
571  xmax = bignum
572  ELSE
573  xmax = xmax*two
574  END IF
575 *
576  IF( notran ) THEN
577 *
578 * Solve A * x = b
579 *
580  DO 120 j = jfirst, jlast, jinc
581 *
582 * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
583 *
584  xj = cabs1( x( j ) )
585  IF( nounit ) THEN
586  tjjs = ab( maind, j )*tscal
587  ELSE
588  tjjs = tscal
589  IF( tscal.EQ.one )
590  $ GO TO 110
591  END IF
592  tjj = cabs1( tjjs )
593  IF( tjj.GT.smlnum ) THEN
594 *
595 * abs(A(j,j)) > SMLNUM:
596 *
597  IF( tjj.LT.one ) THEN
598  IF( xj.GT.tjj*bignum ) THEN
599 *
600 * Scale x by 1/b(j).
601 *
602  rec = one / xj
603  CALL zdscal( n, rec, x, 1 )
604  scale = scale*rec
605  xmax = xmax*rec
606  END IF
607  END IF
608  x( j ) = zladiv( x( j ), tjjs )
609  xj = cabs1( x( j ) )
610  ELSE IF( tjj.GT.zero ) THEN
611 *
612 * 0 < abs(A(j,j)) <= SMLNUM:
613 *
614  IF( xj.GT.tjj*bignum ) THEN
615 *
616 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
617 * to avoid overflow when dividing by A(j,j).
618 *
619  rec = ( tjj*bignum ) / xj
620  IF( cnorm( j ).GT.one ) THEN
621 *
622 * Scale by 1/CNORM(j) to avoid overflow when
623 * multiplying x(j) times column j.
624 *
625  rec = rec / cnorm( j )
626  END IF
627  CALL zdscal( n, rec, x, 1 )
628  scale = scale*rec
629  xmax = xmax*rec
630  END IF
631  x( j ) = zladiv( x( j ), tjjs )
632  xj = cabs1( x( j ) )
633  ELSE
634 *
635 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
636 * scale = 0, and compute a solution to A*x = 0.
637 *
638  DO 100 i = 1, n
639  x( i ) = zero
640  100 CONTINUE
641  x( j ) = one
642  xj = one
643  scale = zero
644  xmax = zero
645  END IF
646  110 CONTINUE
647 *
648 * Scale x if necessary to avoid overflow when adding a
649 * multiple of column j of A.
650 *
651  IF( xj.GT.one ) THEN
652  rec = one / xj
653  IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
654 *
655 * Scale x by 1/(2*abs(x(j))).
656 *
657  rec = rec*half
658  CALL zdscal( n, rec, x, 1 )
659  scale = scale*rec
660  END IF
661  ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
662 *
663 * Scale x by 1/2.
664 *
665  CALL zdscal( n, half, x, 1 )
666  scale = scale*half
667  END IF
668 *
669  IF( upper ) THEN
670  IF( j.GT.1 ) THEN
671 *
672 * Compute the update
673 * x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
674 * x(j)* A(max(1,j-kd):j-1,j)
675 *
676  jlen = min( kd, j-1 )
677  CALL zaxpy( jlen, -x( j )*tscal,
678  $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
679  i = izamax( j-1, x, 1 )
680  xmax = cabs1( x( i ) )
681  END IF
682  ELSE IF( j.LT.n ) THEN
683 *
684 * Compute the update
685 * x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
686 * x(j) * A(j+1:min(j+kd,n),j)
687 *
688  jlen = min( kd, n-j )
689  IF( jlen.GT.0 )
690  $ CALL zaxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
691  $ x( j+1 ), 1 )
692  i = j + izamax( n-j, x( j+1 ), 1 )
693  xmax = cabs1( x( i ) )
694  END IF
695  120 CONTINUE
696 *
697  ELSE IF( lsame( trans, 'T' ) ) THEN
698 *
699 * Solve A**T * x = b
700 *
701  DO 170 j = jfirst, jlast, jinc
702 *
703 * Compute x(j) = b(j) - sum A(k,j)*x(k).
704 * k<>j
705 *
706  xj = cabs1( x( j ) )
707  uscal = tscal
708  rec = one / max( xmax, one )
709  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
710 *
711 * If x(j) could overflow, scale x by 1/(2*XMAX).
712 *
713  rec = rec*half
714  IF( nounit ) THEN
715  tjjs = ab( maind, j )*tscal
716  ELSE
717  tjjs = tscal
718  END IF
719  tjj = cabs1( tjjs )
720  IF( tjj.GT.one ) THEN
721 *
722 * Divide by A(j,j) when scaling x if A(j,j) > 1.
723 *
724  rec = min( one, rec*tjj )
725  uscal = zladiv( uscal, tjjs )
726  END IF
727  IF( rec.LT.one ) THEN
728  CALL zdscal( n, rec, x, 1 )
729  scale = scale*rec
730  xmax = xmax*rec
731  END IF
732  END IF
733 *
734  csumj = zero
735  IF( uscal.EQ.dcmplx( one ) ) THEN
736 *
737 * If the scaling needed for A in the dot product is 1,
738 * call ZDOTU to perform the dot product.
739 *
740  IF( upper ) THEN
741  jlen = min( kd, j-1 )
742  csumj = zdotu( jlen, ab( kd+1-jlen, j ), 1,
743  $ x( j-jlen ), 1 )
744  ELSE
745  jlen = min( kd, n-j )
746  IF( jlen.GT.1 )
747  $ csumj = zdotu( jlen, ab( 2, j ), 1, x( j+1 ),
748  $ 1 )
749  END IF
750  ELSE
751 *
752 * Otherwise, use in-line code for the dot product.
753 *
754  IF( upper ) THEN
755  jlen = min( kd, j-1 )
756  DO 130 i = 1, jlen
757  csumj = csumj + ( ab( kd+i-jlen, j )*uscal )*
758  $ x( j-jlen-1+i )
759  130 CONTINUE
760  ELSE
761  jlen = min( kd, n-j )
762  DO 140 i = 1, jlen
763  csumj = csumj + ( ab( i+1, j )*uscal )*x( j+i )
764  140 CONTINUE
765  END IF
766  END IF
767 *
768  IF( uscal.EQ.dcmplx( tscal ) ) THEN
769 *
770 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
771 * was not used to scale the dotproduct.
772 *
773  x( j ) = x( j ) - csumj
774  xj = cabs1( x( j ) )
775  IF( nounit ) THEN
776 *
777 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
778 *
779  tjjs = ab( maind, j )*tscal
780  ELSE
781  tjjs = tscal
782  IF( tscal.EQ.one )
783  $ GO TO 160
784  END IF
785  tjj = cabs1( tjjs )
786  IF( tjj.GT.smlnum ) THEN
787 *
788 * abs(A(j,j)) > SMLNUM:
789 *
790  IF( tjj.LT.one ) THEN
791  IF( xj.GT.tjj*bignum ) THEN
792 *
793 * Scale X by 1/abs(x(j)).
794 *
795  rec = one / xj
796  CALL zdscal( n, rec, x, 1 )
797  scale = scale*rec
798  xmax = xmax*rec
799  END IF
800  END IF
801  x( j ) = zladiv( x( j ), tjjs )
802  ELSE IF( tjj.GT.zero ) THEN
803 *
804 * 0 < abs(A(j,j)) <= SMLNUM:
805 *
806  IF( xj.GT.tjj*bignum ) THEN
807 *
808 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
809 *
810  rec = ( tjj*bignum ) / xj
811  CALL zdscal( n, rec, x, 1 )
812  scale = scale*rec
813  xmax = xmax*rec
814  END IF
815  x( j ) = zladiv( x( j ), tjjs )
816  ELSE
817 *
818 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
819 * scale = 0 and compute a solution to A**T *x = 0.
820 *
821  DO 150 i = 1, n
822  x( i ) = zero
823  150 CONTINUE
824  x( j ) = one
825  scale = zero
826  xmax = zero
827  END IF
828  160 CONTINUE
829  ELSE
830 *
831 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
832 * product has already been divided by 1/A(j,j).
833 *
834  x( j ) = zladiv( x( j ), tjjs ) - csumj
835  END IF
836  xmax = max( xmax, cabs1( x( j ) ) )
837  170 CONTINUE
838 *
839  ELSE
840 *
841 * Solve A**H * x = b
842 *
843  DO 220 j = jfirst, jlast, jinc
844 *
845 * Compute x(j) = b(j) - sum A(k,j)*x(k).
846 * k<>j
847 *
848  xj = cabs1( x( j ) )
849  uscal = tscal
850  rec = one / max( xmax, one )
851  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
852 *
853 * If x(j) could overflow, scale x by 1/(2*XMAX).
854 *
855  rec = rec*half
856  IF( nounit ) THEN
857  tjjs = dconjg( ab( maind, j ) )*tscal
858  ELSE
859  tjjs = tscal
860  END IF
861  tjj = cabs1( tjjs )
862  IF( tjj.GT.one ) THEN
863 *
864 * Divide by A(j,j) when scaling x if A(j,j) > 1.
865 *
866  rec = min( one, rec*tjj )
867  uscal = zladiv( uscal, tjjs )
868  END IF
869  IF( rec.LT.one ) THEN
870  CALL zdscal( n, rec, x, 1 )
871  scale = scale*rec
872  xmax = xmax*rec
873  END IF
874  END IF
875 *
876  csumj = zero
877  IF( uscal.EQ.dcmplx( one ) ) THEN
878 *
879 * If the scaling needed for A in the dot product is 1,
880 * call ZDOTC to perform the dot product.
881 *
882  IF( upper ) THEN
883  jlen = min( kd, j-1 )
884  csumj = zdotc( jlen, ab( kd+1-jlen, j ), 1,
885  $ x( j-jlen ), 1 )
886  ELSE
887  jlen = min( kd, n-j )
888  IF( jlen.GT.1 )
889  $ csumj = zdotc( jlen, ab( 2, j ), 1, x( j+1 ),
890  $ 1 )
891  END IF
892  ELSE
893 *
894 * Otherwise, use in-line code for the dot product.
895 *
896  IF( upper ) THEN
897  jlen = min( kd, j-1 )
898  DO 180 i = 1, jlen
899  csumj = csumj + ( dconjg( ab( kd+i-jlen, j ) )*
900  $ uscal )*x( j-jlen-1+i )
901  180 CONTINUE
902  ELSE
903  jlen = min( kd, n-j )
904  DO 190 i = 1, jlen
905  csumj = csumj + ( dconjg( ab( i+1, j ) )*uscal )
906  $ *x( j+i )
907  190 CONTINUE
908  END IF
909  END IF
910 *
911  IF( uscal.EQ.dcmplx( tscal ) ) THEN
912 *
913 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
914 * was not used to scale the dotproduct.
915 *
916  x( j ) = x( j ) - csumj
917  xj = cabs1( x( j ) )
918  IF( nounit ) THEN
919 *
920 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
921 *
922  tjjs = dconjg( ab( maind, j ) )*tscal
923  ELSE
924  tjjs = tscal
925  IF( tscal.EQ.one )
926  $ GO TO 210
927  END IF
928  tjj = cabs1( tjjs )
929  IF( tjj.GT.smlnum ) THEN
930 *
931 * abs(A(j,j)) > SMLNUM:
932 *
933  IF( tjj.LT.one ) THEN
934  IF( xj.GT.tjj*bignum ) THEN
935 *
936 * Scale X by 1/abs(x(j)).
937 *
938  rec = one / xj
939  CALL zdscal( n, rec, x, 1 )
940  scale = scale*rec
941  xmax = xmax*rec
942  END IF
943  END IF
944  x( j ) = zladiv( x( j ), tjjs )
945  ELSE IF( tjj.GT.zero ) THEN
946 *
947 * 0 < abs(A(j,j)) <= SMLNUM:
948 *
949  IF( xj.GT.tjj*bignum ) THEN
950 *
951 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
952 *
953  rec = ( tjj*bignum ) / xj
954  CALL zdscal( n, rec, x, 1 )
955  scale = scale*rec
956  xmax = xmax*rec
957  END IF
958  x( j ) = zladiv( x( j ), tjjs )
959  ELSE
960 *
961 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
962 * scale = 0 and compute a solution to A**H *x = 0.
963 *
964  DO 200 i = 1, n
965  x( i ) = zero
966  200 CONTINUE
967  x( j ) = one
968  scale = zero
969  xmax = zero
970  END IF
971  210 CONTINUE
972  ELSE
973 *
974 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
975 * product has already been divided by 1/A(j,j).
976 *
977  x( j ) = zladiv( x( j ), tjjs ) - csumj
978  END IF
979  xmax = max( xmax, cabs1( x( j ) ) )
980  220 CONTINUE
981  END IF
982  scale = scale / tscal
983  END IF
984 *
985 * Scale the column norms by 1/TSCAL for return.
986 *
987  IF( tscal.NE.one ) THEN
988  CALL dscal( n, one / tscal, cnorm, 1 )
989  END IF
990 *
991  RETURN
992 *
993 * End of ZLATBS
994 *
995  END
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
zdscal
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
zaxpy
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
ztbsv
subroutine ztbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
ZTBSV
Definition: ztbsv.f:189
zlatbs
subroutine zlatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
ZLATBS solves a triangular banded system of equations.
Definition: zlatbs.f:243
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79