LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \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)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
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
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
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