LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
slatbs.f
Go to the documentation of this file.
1*> \brief \b SLATBS 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
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slatbs.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slatbs.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slatbs.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLATBS( 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 AB( LDAB, * ), CNORM( * ), X( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SLATBS solves one of the triangular systems
40*>
41*> A *x = s*b or A**T*x = s*b
42*>
43*> with scaling to prevent overflow, where A is an upper or lower
44*> triangular band matrix. Here A**T denotes the transpose of A, x and b
45*> are n-element vectors, and s is a scaling factor, usually less than
46*> or equal to 1, chosen so that the components of x will be less than
47*> the overflow threshold. If the unscaled problem will not cause
48*> overflow, the Level 2 BLAS routine STBSV is called. If the matrix A
49*> is singular (A(j,j) = 0 for some j), then s is set to 0 and a
50*> non-trivial solution to A*x = 0 is returned.
51*> \endverbatim
52*
53* Arguments:
54* ==========
55*
56*> \param[in] UPLO
57*> \verbatim
58*> UPLO is CHARACTER*1
59*> Specifies whether the matrix A is upper or lower triangular.
60*> = 'U': Upper triangular
61*> = 'L': Lower triangular
62*> \endverbatim
63*>
64*> \param[in] TRANS
65*> \verbatim
66*> TRANS is CHARACTER*1
67*> Specifies the operation applied to A.
68*> = 'N': Solve A * x = s*b (No transpose)
69*> = 'T': Solve A**T* x = s*b (Transpose)
70*> = 'C': Solve A**T* x = s*b (Conjugate transpose = Transpose)
71*> \endverbatim
72*>
73*> \param[in] DIAG
74*> \verbatim
75*> DIAG is CHARACTER*1
76*> Specifies whether or not the matrix A is unit triangular.
77*> = 'N': Non-unit triangular
78*> = 'U': Unit triangular
79*> \endverbatim
80*>
81*> \param[in] NORMIN
82*> \verbatim
83*> NORMIN is CHARACTER*1
84*> Specifies whether CNORM has been set or not.
85*> = 'Y': CNORM contains the column norms on entry
86*> = 'N': CNORM is not set on entry. On exit, the norms will
87*> be computed and stored in CNORM.
88*> \endverbatim
89*>
90*> \param[in] N
91*> \verbatim
92*> N is INTEGER
93*> The order of the matrix A. N >= 0.
94*> \endverbatim
95*>
96*> \param[in] KD
97*> \verbatim
98*> KD is INTEGER
99*> The number of subdiagonals or superdiagonals in the
100*> triangular matrix A. KD >= 0.
101*> \endverbatim
102*>
103*> \param[in] AB
104*> \verbatim
105*> AB is REAL array, dimension (LDAB,N)
106*> The upper or lower triangular band matrix A, stored in the
107*> first KD+1 rows of the array. The j-th column of A is stored
108*> in the j-th column of the array AB as follows:
109*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
110*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
111*> \endverbatim
112*>
113*> \param[in] LDAB
114*> \verbatim
115*> LDAB is INTEGER
116*> The leading dimension of the array AB. LDAB >= KD+1.
117*> \endverbatim
118*>
119*> \param[in,out] X
120*> \verbatim
121*> X is REAL array, dimension (N)
122*> On entry, the right hand side b of the triangular system.
123*> On exit, X is overwritten by the solution vector x.
124*> \endverbatim
125*>
126*> \param[out] SCALE
127*> \verbatim
128*> SCALE is REAL
129*> The scaling factor s for the triangular system
130*> A * x = s*b or A**T* x = s*b.
131*> If SCALE = 0, the matrix A is singular or badly scaled, and
132*> the vector x is an exact or approximate solution to A*x = 0.
133*> \endverbatim
134*>
135*> \param[in,out] CNORM
136*> \verbatim
137*> CNORM is REAL array, dimension (N)
138*>
139*> If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
140*> contains the norm of the off-diagonal part of the j-th column
141*> of A. If TRANS = 'N', CNORM(j) must be greater than or equal
142*> to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
143*> must be greater than or equal to the 1-norm.
144*>
145*> If NORMIN = 'N', CNORM is an output argument and CNORM(j)
146*> returns the 1-norm of the offdiagonal part of the j-th column
147*> of A.
148*> \endverbatim
149*>
150*> \param[out] INFO
151*> \verbatim
152*> INFO is INTEGER
153*> = 0: successful exit
154*> < 0: if INFO = -k, the k-th argument had an illegal value
155*> \endverbatim
156*
157* Authors:
158* ========
159*
160*> \author Univ. of Tennessee
161*> \author Univ. of California Berkeley
162*> \author Univ. of Colorado Denver
163*> \author NAG Ltd.
164*
165*> \ingroup realOTHERauxiliary
166*
167*> \par Further Details:
168* =====================
169*>
170*> \verbatim
171*>
172*> A rough bound on x is computed; if that is less than overflow, STBSV
173*> is called, otherwise, specific code is used which checks for possible
174*> overflow or divide-by-zero at every operation.
175*>
176*> A columnwise scheme is used for solving A*x = b. The basic algorithm
177*> if A is lower triangular is
178*>
179*> x[1:n] := b[1:n]
180*> for j = 1, ..., n
181*> x(j) := x(j) / A(j,j)
182*> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
183*> end
184*>
185*> Define bounds on the components of x after j iterations of the loop:
186*> M(j) = bound on x[1:j]
187*> G(j) = bound on x[j+1:n]
188*> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
189*>
190*> Then for iteration j+1 we have
191*> M(j+1) <= G(j) / | A(j+1,j+1) |
192*> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
193*> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
194*>
195*> where CNORM(j+1) is greater than or equal to the infinity-norm of
196*> column j+1 of A, not counting the diagonal. Hence
197*>
198*> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
199*> 1<=i<=j
200*> and
201*>
202*> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
203*> 1<=i< j
204*>
205*> Since |x(j)| <= M(j), we use the Level 2 BLAS routine STBSV if the
206*> reciprocal of the largest M(j), j=1,..,n, is larger than
207*> max(underflow, 1/overflow).
208*>
209*> The bound on x(j) is also used to determine when a step in the
210*> columnwise method can be performed without fear of overflow. If
211*> the computed bound is greater than a large constant, x is scaled to
212*> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
213*> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
214*>
215*> Similarly, a row-wise scheme is used to solve A**T*x = b. The basic
216*> algorithm for A upper triangular is
217*>
218*> for j = 1, ..., n
219*> x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
220*> end
221*>
222*> We simultaneously compute two bounds
223*> G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
224*> M(j) = bound on x(i), 1<=i<=j
225*>
226*> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
227*> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
228*> Then the bound on x(j) is
229*>
230*> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
231*>
232*> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
233*> 1<=i<=j
234*>
235*> and we can safely call STBSV if 1/M(n) and 1/G(n) are both greater
236*> than max(underflow, 1/overflow).
237*> \endverbatim
238*>
239* =====================================================================
240 SUBROUTINE slatbs( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
241 \$ SCALE, CNORM, INFO )
242*
243* -- LAPACK auxiliary routine --
244* -- LAPACK is a software package provided by Univ. of Tennessee, --
245* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
246*
247* .. Scalar Arguments ..
248 CHARACTER DIAG, NORMIN, TRANS, UPLO
249 INTEGER INFO, KD, LDAB, N
250 REAL SCALE
251* ..
252* .. Array Arguments ..
253 REAL AB( LDAB, * ), CNORM( * ), X( * )
254* ..
255*
256* =====================================================================
257*
258* .. Parameters ..
259 REAL ZERO, HALF, ONE
260 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0 )
261* ..
262* .. Local Scalars ..
263 LOGICAL NOTRAN, NOUNIT, UPPER
264 INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
265 REAL BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
266 \$ tmax, tscal, uscal, xbnd, xj, xmax
267* ..
268* .. External Functions ..
269 LOGICAL LSAME
270 INTEGER ISAMAX
271 REAL SASUM, SDOT, SLAMCH
272 EXTERNAL lsame, isamax, sasum, sdot, slamch
273* ..
274* .. External Subroutines ..
275 EXTERNAL saxpy, sscal, stbsv, xerbla
276* ..
277* .. Intrinsic Functions ..
278 INTRINSIC abs, max, min
279* ..
280* .. Executable Statements ..
281*
282 info = 0
283 upper = lsame( uplo, 'U' )
284 notran = lsame( trans, 'N' )
285 nounit = lsame( diag, 'N' )
286*
287* Test the input parameters.
288*
289 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
290 info = -1
291 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
292 \$ lsame( trans, 'C' ) ) THEN
293 info = -2
294 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
295 info = -3
296 ELSE IF( .NOT.lsame( normin, 'Y' ) .AND. .NOT.
297 \$ lsame( normin, 'N' ) ) THEN
298 info = -4
299 ELSE IF( n.LT.0 ) THEN
300 info = -5
301 ELSE IF( kd.LT.0 ) THEN
302 info = -6
303 ELSE IF( ldab.LT.kd+1 ) THEN
304 info = -8
305 END IF
306 IF( info.NE.0 ) THEN
307 CALL xerbla( 'SLATBS', -info )
308 RETURN
309 END IF
310*
311* Quick return if possible
312*
313 scale = one
314 IF( n.EQ.0 )
315 \$ RETURN
316*
317* Determine machine dependent parameters to control overflow.
318*
319 smlnum = slamch( 'Safe minimum' ) / slamch( 'Precision' )
320 bignum = one / smlnum
321*
322 IF( lsame( normin, 'N' ) ) THEN
323*
324* Compute the 1-norm of each column, not including the diagonal.
325*
326 IF( upper ) THEN
327*
328* A is upper triangular.
329*
330 DO 10 j = 1, n
331 jlen = min( kd, j-1 )
332 cnorm( j ) = sasum( jlen, ab( kd+1-jlen, j ), 1 )
333 10 CONTINUE
334 ELSE
335*
336* A is lower triangular.
337*
338 DO 20 j = 1, n
339 jlen = min( kd, n-j )
340 IF( jlen.GT.0 ) THEN
341 cnorm( j ) = sasum( jlen, ab( 2, j ), 1 )
342 ELSE
343 cnorm( j ) = zero
344 END IF
345 20 CONTINUE
346 END IF
347 END IF
348*
349* Scale the column norms by TSCAL if the maximum element in CNORM is
350* greater than BIGNUM.
351*
352 imax = isamax( n, cnorm, 1 )
353 tmax = cnorm( imax )
354 IF( tmax.LE.bignum ) THEN
355 tscal = one
356 ELSE
357 tscal = one / ( smlnum*tmax )
358 CALL sscal( n, tscal, cnorm, 1 )
359 END IF
360*
361* Compute a bound on the computed solution vector to see if the
362* Level 2 BLAS routine STBSV can be used.
363*
364 j = isamax( n, x, 1 )
365 xmax = abs( x( j ) )
366 xbnd = xmax
367 IF( notran ) THEN
368*
369* Compute the growth in A * x = b.
370*
371 IF( upper ) THEN
372 jfirst = n
373 jlast = 1
374 jinc = -1
375 maind = kd + 1
376 ELSE
377 jfirst = 1
378 jlast = n
379 jinc = 1
380 maind = 1
381 END IF
382*
383 IF( tscal.NE.one ) THEN
384 grow = zero
385 GO TO 50
386 END IF
387*
388 IF( nounit ) THEN
389*
390* A is non-unit triangular.
391*
392* Compute GROW = 1/G(j) and XBND = 1/M(j).
393* Initially, G(0) = max{x(i), i=1,...,n}.
394*
395 grow = one / max( xbnd, smlnum )
396 xbnd = grow
397 DO 30 j = jfirst, jlast, jinc
398*
399* Exit the loop if the growth factor is too small.
400*
401 IF( grow.LE.smlnum )
402 \$ GO TO 50
403*
404* M(j) = G(j-1) / abs(A(j,j))
405*
406 tjj = abs( ab( maind, j ) )
407 xbnd = min( xbnd, min( one, tjj )*grow )
408 IF( tjj+cnorm( j ).GE.smlnum ) THEN
409*
410* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
411*
412 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
413 ELSE
414*
415* G(j) could overflow, set GROW to 0.
416*
417 grow = zero
418 END IF
419 30 CONTINUE
420 grow = xbnd
421 ELSE
422*
423* A is unit triangular.
424*
425* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
426*
427 grow = min( one, one / max( xbnd, smlnum ) )
428 DO 40 j = jfirst, jlast, jinc
429*
430* Exit the loop if the growth factor is too small.
431*
432 IF( grow.LE.smlnum )
433 \$ GO TO 50
434*
435* G(j) = G(j-1)*( 1 + CNORM(j) )
436*
437 grow = grow*( one / ( one+cnorm( j ) ) )
438 40 CONTINUE
439 END IF
440 50 CONTINUE
441*
442 ELSE
443*
444* Compute the growth in A**T * x = b.
445*
446 IF( upper ) THEN
447 jfirst = 1
448 jlast = n
449 jinc = 1
450 maind = kd + 1
451 ELSE
452 jfirst = n
453 jlast = 1
454 jinc = -1
455 maind = 1
456 END IF
457*
458 IF( tscal.NE.one ) THEN
459 grow = zero
460 GO TO 80
461 END IF
462*
463 IF( nounit ) THEN
464*
465* A is non-unit triangular.
466*
467* Compute GROW = 1/G(j) and XBND = 1/M(j).
468* Initially, M(0) = max{x(i), i=1,...,n}.
469*
470 grow = one / max( xbnd, smlnum )
471 xbnd = grow
472 DO 60 j = jfirst, jlast, jinc
473*
474* Exit the loop if the growth factor is too small.
475*
476 IF( grow.LE.smlnum )
477 \$ GO TO 80
478*
479* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
480*
481 xj = one + cnorm( j )
482 grow = min( grow, xbnd / xj )
483*
484* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
485*
486 tjj = abs( ab( maind, j ) )
487 IF( xj.GT.tjj )
488 \$ xbnd = xbnd*( tjj / xj )
489 60 CONTINUE
490 grow = min( grow, xbnd )
491 ELSE
492*
493* A is unit triangular.
494*
495* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
496*
497 grow = min( one, one / max( xbnd, smlnum ) )
498 DO 70 j = jfirst, jlast, jinc
499*
500* Exit the loop if the growth factor is too small.
501*
502 IF( grow.LE.smlnum )
503 \$ GO TO 80
504*
505* G(j) = ( 1 + CNORM(j) )*G(j-1)
506*
507 xj = one + cnorm( j )
508 grow = grow / xj
509 70 CONTINUE
510 END IF
511 80 CONTINUE
512 END IF
513*
514 IF( ( grow*tscal ).GT.smlnum ) THEN
515*
516* Use the Level 2 BLAS solve if the reciprocal of the bound on
517* elements of X is not too small.
518*
519 CALL stbsv( uplo, trans, diag, n, kd, ab, ldab, x, 1 )
520 ELSE
521*
522* Use a Level 1 BLAS solve, scaling intermediate results.
523*
524 IF( xmax.GT.bignum ) THEN
525*
526* Scale X so that its components are less than or equal to
527* BIGNUM in absolute value.
528*
529 scale = bignum / xmax
530 CALL sscal( n, scale, x, 1 )
531 xmax = bignum
532 END IF
533*
534 IF( notran ) THEN
535*
536* Solve A * x = b
537*
538 DO 100 j = jfirst, jlast, jinc
539*
540* Compute x(j) = b(j) / A(j,j), scaling x if necessary.
541*
542 xj = abs( x( j ) )
543 IF( nounit ) THEN
544 tjjs = ab( maind, j )*tscal
545 ELSE
546 tjjs = tscal
547 IF( tscal.EQ.one )
548 \$ GO TO 95
549 END IF
550 tjj = abs( tjjs )
551 IF( tjj.GT.smlnum ) THEN
552*
553* abs(A(j,j)) > SMLNUM:
554*
555 IF( tjj.LT.one ) THEN
556 IF( xj.GT.tjj*bignum ) THEN
557*
558* Scale x by 1/b(j).
559*
560 rec = one / xj
561 CALL sscal( n, rec, x, 1 )
562 scale = scale*rec
563 xmax = xmax*rec
564 END IF
565 END IF
566 x( j ) = x( j ) / tjjs
567 xj = abs( x( j ) )
568 ELSE IF( tjj.GT.zero ) THEN
569*
570* 0 < abs(A(j,j)) <= SMLNUM:
571*
572 IF( xj.GT.tjj*bignum ) THEN
573*
574* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
575* to avoid overflow when dividing by A(j,j).
576*
577 rec = ( tjj*bignum ) / xj
578 IF( cnorm( j ).GT.one ) THEN
579*
580* Scale by 1/CNORM(j) to avoid overflow when
581* multiplying x(j) times column j.
582*
583 rec = rec / cnorm( j )
584 END IF
585 CALL sscal( n, rec, x, 1 )
586 scale = scale*rec
587 xmax = xmax*rec
588 END IF
589 x( j ) = x( j ) / tjjs
590 xj = abs( x( j ) )
591 ELSE
592*
593* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
594* scale = 0, and compute a solution to A*x = 0.
595*
596 DO 90 i = 1, n
597 x( i ) = zero
598 90 CONTINUE
599 x( j ) = one
600 xj = one
601 scale = zero
602 xmax = zero
603 END IF
604 95 CONTINUE
605*
606* Scale x if necessary to avoid overflow when adding a
607* multiple of column j of A.
608*
609 IF( xj.GT.one ) THEN
610 rec = one / xj
611 IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
612*
613* Scale x by 1/(2*abs(x(j))).
614*
615 rec = rec*half
616 CALL sscal( n, rec, x, 1 )
617 scale = scale*rec
618 END IF
619 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
620*
621* Scale x by 1/2.
622*
623 CALL sscal( n, half, x, 1 )
624 scale = scale*half
625 END IF
626*
627 IF( upper ) THEN
628 IF( j.GT.1 ) THEN
629*
630* Compute the update
631* x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
632* x(j)* A(max(1,j-kd):j-1,j)
633*
634 jlen = min( kd, j-1 )
635 CALL saxpy( jlen, -x( j )*tscal,
636 \$ ab( kd+1-jlen, j ), 1, x( j-jlen ), 1 )
637 i = isamax( j-1, x, 1 )
638 xmax = abs( x( i ) )
639 END IF
640 ELSE IF( j.LT.n ) THEN
641*
642* Compute the update
643* x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
644* x(j) * A(j+1:min(j+kd,n),j)
645*
646 jlen = min( kd, n-j )
647 IF( jlen.GT.0 )
648 \$ CALL saxpy( jlen, -x( j )*tscal, ab( 2, j ), 1,
649 \$ x( j+1 ), 1 )
650 i = j + isamax( n-j, x( j+1 ), 1 )
651 xmax = abs( x( i ) )
652 END IF
653 100 CONTINUE
654*
655 ELSE
656*
657* Solve A**T * x = b
658*
659 DO 140 j = jfirst, jlast, jinc
660*
661* Compute x(j) = b(j) - sum A(k,j)*x(k).
662* k<>j
663*
664 xj = abs( x( j ) )
665 uscal = tscal
666 rec = one / max( xmax, one )
667 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
668*
669* If x(j) could overflow, scale x by 1/(2*XMAX).
670*
671 rec = rec*half
672 IF( nounit ) THEN
673 tjjs = ab( maind, j )*tscal
674 ELSE
675 tjjs = tscal
676 END IF
677 tjj = abs( tjjs )
678 IF( tjj.GT.one ) THEN
679*
680* Divide by A(j,j) when scaling x if A(j,j) > 1.
681*
682 rec = min( one, rec*tjj )
683 uscal = uscal / tjjs
684 END IF
685 IF( rec.LT.one ) THEN
686 CALL sscal( n, rec, x, 1 )
687 scale = scale*rec
688 xmax = xmax*rec
689 END IF
690 END IF
691*
692 sumj = zero
693 IF( uscal.EQ.one ) THEN
694*
695* If the scaling needed for A in the dot product is 1,
696* call SDOT to perform the dot product.
697*
698 IF( upper ) THEN
699 jlen = min( kd, j-1 )
700 sumj = sdot( jlen, ab( kd+1-jlen, j ), 1,
701 \$ x( j-jlen ), 1 )
702 ELSE
703 jlen = min( kd, n-j )
704 IF( jlen.GT.0 )
705 \$ sumj = sdot( jlen, ab( 2, j ), 1, x( j+1 ), 1 )
706 END IF
707 ELSE
708*
709* Otherwise, use in-line code for the dot product.
710*
711 IF( upper ) THEN
712 jlen = min( kd, j-1 )
713 DO 110 i = 1, jlen
714 sumj = sumj + ( ab( kd+i-jlen, j )*uscal )*
715 \$ x( j-jlen-1+i )
716 110 CONTINUE
717 ELSE
718 jlen = min( kd, n-j )
719 DO 120 i = 1, jlen
720 sumj = sumj + ( ab( i+1, j )*uscal )*x( j+i )
721 120 CONTINUE
722 END IF
723 END IF
724*
725 IF( uscal.EQ.tscal ) THEN
726*
727* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
728* was not used to scale the dotproduct.
729*
730 x( j ) = x( j ) - sumj
731 xj = abs( x( j ) )
732 IF( nounit ) THEN
733*
734* Compute x(j) = x(j) / A(j,j), scaling if necessary.
735*
736 tjjs = ab( maind, j )*tscal
737 ELSE
738 tjjs = tscal
739 IF( tscal.EQ.one )
740 \$ GO TO 135
741 END IF
742 tjj = abs( tjjs )
743 IF( tjj.GT.smlnum ) THEN
744*
745* abs(A(j,j)) > SMLNUM:
746*
747 IF( tjj.LT.one ) THEN
748 IF( xj.GT.tjj*bignum ) THEN
749*
750* Scale X by 1/abs(x(j)).
751*
752 rec = one / xj
753 CALL sscal( n, rec, x, 1 )
754 scale = scale*rec
755 xmax = xmax*rec
756 END IF
757 END IF
758 x( j ) = x( j ) / tjjs
759 ELSE IF( tjj.GT.zero ) THEN
760*
761* 0 < abs(A(j,j)) <= SMLNUM:
762*
763 IF( xj.GT.tjj*bignum ) THEN
764*
765* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
766*
767 rec = ( tjj*bignum ) / xj
768 CALL sscal( n, rec, x, 1 )
769 scale = scale*rec
770 xmax = xmax*rec
771 END IF
772 x( j ) = x( j ) / tjjs
773 ELSE
774*
775* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
776* scale = 0, and compute a solution to A**T*x = 0.
777*
778 DO 130 i = 1, n
779 x( i ) = zero
780 130 CONTINUE
781 x( j ) = one
782 scale = zero
783 xmax = zero
784 END IF
785 135 CONTINUE
786 ELSE
787*
788* Compute x(j) := x(j) / A(j,j) - sumj if the dot
789* product has already been divided by 1/A(j,j).
790*
791 x( j ) = x( j ) / tjjs - sumj
792 END IF
793 xmax = max( xmax, abs( x( j ) ) )
794 140 CONTINUE
795 END IF
796 scale = scale / tscal
797 END IF
798*
799* Scale the column norms by 1/TSCAL for return.
800*
801 IF( tscal.NE.one ) THEN
802 CALL sscal( n, one / tscal, cnorm, 1 )
803 END IF
804*
805 RETURN
806*
807* End of SLATBS
808*
809 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
SLATBS solves a triangular banded system of equations.
Definition: slatbs.f:242
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine stbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
STBSV
Definition: stbsv.f:189