OpenRadioss 2025.1.11
OpenRadioss project
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 complexOTHERauxiliary
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, slabad, 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 IF( n.EQ.0 )
328 $ RETURN
329*
330* Determine machine dependent parameters to control overflow.
331*
332 smlnum = slamch( 'Safe minimum' )
333 bignum = one / smlnum
334 CALL slabad( smlnum, bignum )
335 smlnum = smlnum / slamch( '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 ) = scasum( 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 ) = scasum( 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 = isamax( 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 sscal( 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 CTBSV 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 ctbsv( 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 csscal( 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 110 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 105
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 csscal( n, rec, x, 1 )
604 scale = scale*rec
605 xmax = xmax*rec
606 END IF
607 END IF
608 x( j ) = cladiv( 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 csscal( n, rec, x, 1 )
628 scale = scale*rec
629 xmax = xmax*rec
630 END IF
631 x( j ) = cladiv( 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 105 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 csscal( 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 csscal( 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 caxpy( jlen, -x( j )*tscal,
678 $ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
679 i = icamax( 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 caxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
691 $ x( j+1 ), 1 )
692 i = j + icamax( n-j, x( j+1 ), 1 )
693 xmax = cabs1( x( i ) )
694 END IF
695 110 CONTINUE
696*
697 ELSE IF( lsame( trans, 'T' ) ) THEN
698*
699* Solve A**T * x = b
700*
701 DO 150 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 = cladiv( uscal, tjjs )
726 END IF
727 IF( rec.LT.one ) THEN
728 CALL csscal( 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.cmplx( one ) ) THEN
736*
737* If the scaling needed for A in the dot product is 1,
738* call CDOTU to perform the dot product.
739*
740 IF( upper ) THEN
741 jlen = min( kd, j-1 )
742 csumj = cdotu( 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 = cdotu( 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 120 i = 1, jlen
757 csumj = csumj + ( ab( kd+i-jlen, j )*uscal )*
758 $ x( j-jlen-1+i )
759 120 CONTINUE
760 ELSE
761 jlen = min( kd, n-j )
762 DO 130 i = 1, jlen
763 csumj = csumj + ( ab( i+1, j )*uscal )*x( j+i )
764 130 CONTINUE
765 END IF
766 END IF
767*
768 IF( uscal.EQ.cmplx( 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 145
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 csscal( n, rec, x, 1 )
797 scale = scale*rec
798 xmax = xmax*rec
799 END IF
800 END IF
801 x( j ) = cladiv( 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 csscal( n, rec, x, 1 )
812 scale = scale*rec
813 xmax = xmax*rec
814 END IF
815 x( j ) = cladiv( 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 140 i = 1, n
822 x( i ) = zero
823 140 CONTINUE
824 x( j ) = one
825 scale = zero
826 xmax = zero
827 END IF
828 145 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 ) = cladiv( x( j ), tjjs ) - csumj
835 END IF
836 xmax = max( xmax, cabs1( x( j ) ) )
837 150 CONTINUE
838*
839 ELSE
840*
841* Solve A**H * x = b
842*
843 DO 190 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 = conjg( 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 = cladiv( uscal, tjjs )
868 END IF
869 IF( rec.LT.one ) THEN
870 CALL csscal( 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.cmplx( one ) ) THEN
878*
879* If the scaling needed for A in the dot product is 1,
880* call CDOTC to perform the dot product.
881*
882 IF( upper ) THEN
883 jlen = min( kd, j-1 )
884 csumj = cdotc( 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 = cdotc( 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 160 i = 1, jlen
899 csumj = csumj + ( conjg( ab( kd+i-jlen, j ) )*
900 $ uscal )*x( j-jlen-1+i )
901 160 CONTINUE
902 ELSE
903 jlen = min( kd, n-j )
904 DO 170 i = 1, jlen
905 csumj = csumj + ( conjg( ab( i+1, j ) )*uscal )*
906 $ x( j+i )
907 170 CONTINUE
908 END IF
909 END IF
910*
911 IF( uscal.EQ.cmplx( 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 = conjg( ab( maind, j ) )*tscal
923 ELSE
924 tjjs = tscal
925 IF( tscal.EQ.one )
926 $ GO TO 185
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 csscal( n, rec, x, 1 )
940 scale = scale*rec
941 xmax = xmax*rec
942 END IF
943 END IF
944 x( j ) = cladiv( 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 csscal( n, rec, x, 1 )
955 scale = scale*rec
956 xmax = xmax*rec
957 END IF
958 x( j ) = cladiv( 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 180 i = 1, n
965 x( i ) = zero
966 180 CONTINUE
967 x( j ) = one
968 scale = zero
969 xmax = zero
970 END IF
971 185 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 ) = cladiv( x( j ), tjjs ) - csumj
978 END IF
979 xmax = max( xmax, cabs1( x( j ) ) )
980 190 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 sscal( n, one / tscal, cnorm, 1 )
989 END IF
990*
991 RETURN
992*
993* End of CLATBS
994*
995 END
cmplx
float cmplx[2]
Definition pblas.h:136
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
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
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
ctbsv
subroutine ctbsv(uplo, trans, diag, n, k, a, lda, x, incx)
CTBSV
Definition ctbsv.f:189
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21