![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlatbs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dlatbs.f CVS log | Up to [local] / rpl / lapack / lapack |
1.9 bertrand 1: *> \brief \b DLATBS
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLATBS + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlatbs.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlatbs.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlatbs.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLATBS( 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 AB( LDAB, * ), CNORM( * ), X( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DLATBS 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 DTBSV 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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
1.11 bertrand 137: *> CNORM is DOUBLE PRECISION array, dimension (N)
1.9 bertrand 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: *
1.11 bertrand 165: *> \date April 2012
1.9 bertrand 166: *
167: *> \ingroup doubleOTHERauxiliary
168: *
169: *> \par Further Details:
170: * =====================
171: *>
172: *> \verbatim
173: *>
174: *> A rough bound on x is computed; if that is less than overflow, DTBSV
175: *> is called, otherwise, specific code is used which checks for possible
176: *> overflow or divide-by-zero at every operation.
177: *>
178: *> A columnwise scheme is used for solving A*x = b. The basic algorithm
179: *> if A is lower triangular is
180: *>
181: *> x[1:n] := b[1:n]
182: *> for j = 1, ..., n
183: *> x(j) := x(j) / A(j,j)
184: *> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
185: *> end
186: *>
187: *> Define bounds on the components of x after j iterations of the loop:
188: *> M(j) = bound on x[1:j]
189: *> G(j) = bound on x[j+1:n]
190: *> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
191: *>
192: *> Then for iteration j+1 we have
193: *> M(j+1) <= G(j) / | A(j+1,j+1) |
194: *> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
195: *> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
196: *>
197: *> where CNORM(j+1) is greater than or equal to the infinity-norm of
198: *> column j+1 of A, not counting the diagonal. Hence
199: *>
200: *> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
201: *> 1<=i<=j
202: *> and
203: *>
204: *> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
205: *> 1<=i< j
206: *>
207: *> Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTBSV if the
208: *> reciprocal of the largest M(j), j=1,..,n, is larger than
209: *> max(underflow, 1/overflow).
210: *>
211: *> The bound on x(j) is also used to determine when a step in the
212: *> columnwise method can be performed without fear of overflow. If
213: *> the computed bound is greater than a large constant, x is scaled to
214: *> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
215: *> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
216: *>
217: *> Similarly, a row-wise scheme is used to solve A**T*x = b. The basic
218: *> algorithm for A upper triangular is
219: *>
220: *> for j = 1, ..., n
221: *> x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
222: *> end
223: *>
224: *> We simultaneously compute two bounds
225: *> G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
226: *> M(j) = bound on x(i), 1<=i<=j
227: *>
228: *> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
229: *> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
230: *> Then the bound on x(j) is
231: *>
232: *> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
233: *>
234: *> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
235: *> 1<=i<=j
236: *>
237: *> and we can safely call DTBSV if 1/M(n) and 1/G(n) are both greater
238: *> than max(underflow, 1/overflow).
239: *> \endverbatim
240: *>
241: * =====================================================================
1.1 bertrand 242: SUBROUTINE DLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
243: $ SCALE, CNORM, INFO )
244: *
1.11 bertrand 245: * -- LAPACK auxiliary routine (version 3.4.1) --
1.1 bertrand 246: * -- LAPACK is a software package provided by Univ. of Tennessee, --
247: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11 bertrand 248: * April 2012
1.1 bertrand 249: *
250: * .. Scalar Arguments ..
251: CHARACTER DIAG, NORMIN, TRANS, UPLO
252: INTEGER INFO, KD, LDAB, N
253: DOUBLE PRECISION SCALE
254: * ..
255: * .. Array Arguments ..
256: DOUBLE PRECISION AB( LDAB, * ), CNORM( * ), X( * )
257: * ..
258: *
259: * =====================================================================
260: *
261: * .. Parameters ..
262: DOUBLE PRECISION ZERO, HALF, ONE
263: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.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, SUMJ, TJJ, TJJS,
269: $ TMAX, TSCAL, USCAL, XBND, XJ, XMAX
270: * ..
271: * .. External Functions ..
272: LOGICAL LSAME
273: INTEGER IDAMAX
274: DOUBLE PRECISION DASUM, DDOT, DLAMCH
275: EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
276: * ..
277: * .. External Subroutines ..
278: EXTERNAL DAXPY, DSCAL, DTBSV, XERBLA
279: * ..
280: * .. Intrinsic Functions ..
281: INTRINSIC ABS, MAX, MIN
282: * ..
283: * .. Executable Statements ..
284: *
285: INFO = 0
286: UPPER = LSAME( UPLO, 'U' )
287: NOTRAN = LSAME( TRANS, 'N' )
288: NOUNIT = LSAME( DIAG, 'N' )
289: *
290: * Test the input parameters.
291: *
292: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
293: INFO = -1
294: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
295: $ LSAME( TRANS, 'C' ) ) THEN
296: INFO = -2
297: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
298: INFO = -3
299: ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
300: $ LSAME( NORMIN, 'N' ) ) THEN
301: INFO = -4
302: ELSE IF( N.LT.0 ) THEN
303: INFO = -5
304: ELSE IF( KD.LT.0 ) THEN
305: INFO = -6
306: ELSE IF( LDAB.LT.KD+1 ) THEN
307: INFO = -8
308: END IF
309: IF( INFO.NE.0 ) THEN
310: CALL XERBLA( 'DLATBS', -INFO )
311: RETURN
312: END IF
313: *
314: * Quick return if possible
315: *
316: IF( N.EQ.0 )
317: $ RETURN
318: *
319: * Determine machine dependent parameters to control overflow.
320: *
321: SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
322: BIGNUM = ONE / SMLNUM
323: SCALE = ONE
324: *
325: IF( LSAME( NORMIN, 'N' ) ) THEN
326: *
327: * Compute the 1-norm of each column, not including the diagonal.
328: *
329: IF( UPPER ) THEN
330: *
331: * A is upper triangular.
332: *
333: DO 10 J = 1, N
334: JLEN = MIN( KD, J-1 )
335: CNORM( J ) = DASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
336: 10 CONTINUE
337: ELSE
338: *
339: * A is lower triangular.
340: *
341: DO 20 J = 1, N
342: JLEN = MIN( KD, N-J )
343: IF( JLEN.GT.0 ) THEN
344: CNORM( J ) = DASUM( JLEN, AB( 2, J ), 1 )
345: ELSE
346: CNORM( J ) = ZERO
347: END IF
348: 20 CONTINUE
349: END IF
350: END IF
351: *
352: * Scale the column norms by TSCAL if the maximum element in CNORM is
353: * greater than BIGNUM.
354: *
355: IMAX = IDAMAX( N, CNORM, 1 )
356: TMAX = CNORM( IMAX )
357: IF( TMAX.LE.BIGNUM ) THEN
358: TSCAL = ONE
359: ELSE
360: TSCAL = ONE / ( SMLNUM*TMAX )
361: CALL DSCAL( N, TSCAL, CNORM, 1 )
362: END IF
363: *
364: * Compute a bound on the computed solution vector to see if the
365: * Level 2 BLAS routine DTBSV can be used.
366: *
367: J = IDAMAX( N, X, 1 )
368: XMAX = ABS( X( J ) )
369: XBND = XMAX
370: IF( NOTRAN ) THEN
371: *
372: * Compute the growth in A * x = b.
373: *
374: IF( UPPER ) THEN
375: JFIRST = N
376: JLAST = 1
377: JINC = -1
378: MAIND = KD + 1
379: ELSE
380: JFIRST = 1
381: JLAST = N
382: JINC = 1
383: MAIND = 1
384: END IF
385: *
386: IF( TSCAL.NE.ONE ) THEN
387: GROW = ZERO
388: GO TO 50
389: END IF
390: *
391: IF( NOUNIT ) THEN
392: *
393: * A is non-unit triangular.
394: *
395: * Compute GROW = 1/G(j) and XBND = 1/M(j).
396: * Initially, G(0) = max{x(i), i=1,...,n}.
397: *
398: GROW = ONE / MAX( XBND, SMLNUM )
399: XBND = GROW
400: DO 30 J = JFIRST, JLAST, JINC
401: *
402: * Exit the loop if the growth factor is too small.
403: *
404: IF( GROW.LE.SMLNUM )
405: $ GO TO 50
406: *
407: * M(j) = G(j-1) / abs(A(j,j))
408: *
409: TJJ = ABS( AB( MAIND, J ) )
410: XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
411: IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
412: *
413: * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
414: *
415: GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
416: ELSE
417: *
418: * G(j) could overflow, set GROW to 0.
419: *
420: GROW = ZERO
421: END IF
422: 30 CONTINUE
423: GROW = XBND
424: ELSE
425: *
426: * A is unit triangular.
427: *
428: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
429: *
430: GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
431: DO 40 J = JFIRST, JLAST, JINC
432: *
433: * Exit the loop if the growth factor is too small.
434: *
435: IF( GROW.LE.SMLNUM )
436: $ GO TO 50
437: *
438: * G(j) = G(j-1)*( 1 + CNORM(j) )
439: *
440: GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
441: 40 CONTINUE
442: END IF
443: 50 CONTINUE
444: *
445: ELSE
446: *
1.8 bertrand 447: * Compute the growth in A**T * x = b.
1.1 bertrand 448: *
449: IF( UPPER ) THEN
450: JFIRST = 1
451: JLAST = N
452: JINC = 1
453: MAIND = KD + 1
454: ELSE
455: JFIRST = N
456: JLAST = 1
457: JINC = -1
458: MAIND = 1
459: END IF
460: *
461: IF( TSCAL.NE.ONE ) THEN
462: GROW = ZERO
463: GO TO 80
464: END IF
465: *
466: IF( NOUNIT ) THEN
467: *
468: * A is non-unit triangular.
469: *
470: * Compute GROW = 1/G(j) and XBND = 1/M(j).
471: * Initially, M(0) = max{x(i), i=1,...,n}.
472: *
473: GROW = ONE / MAX( XBND, SMLNUM )
474: XBND = GROW
475: DO 60 J = JFIRST, JLAST, JINC
476: *
477: * Exit the loop if the growth factor is too small.
478: *
479: IF( GROW.LE.SMLNUM )
480: $ GO TO 80
481: *
482: * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
483: *
484: XJ = ONE + CNORM( J )
485: GROW = MIN( GROW, XBND / XJ )
486: *
487: * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
488: *
489: TJJ = ABS( AB( MAIND, J ) )
490: IF( XJ.GT.TJJ )
491: $ XBND = XBND*( TJJ / XJ )
492: 60 CONTINUE
493: GROW = MIN( GROW, XBND )
494: ELSE
495: *
496: * A is unit triangular.
497: *
498: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
499: *
500: GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
501: DO 70 J = JFIRST, JLAST, JINC
502: *
503: * Exit the loop if the growth factor is too small.
504: *
505: IF( GROW.LE.SMLNUM )
506: $ GO TO 80
507: *
508: * G(j) = ( 1 + CNORM(j) )*G(j-1)
509: *
510: XJ = ONE + CNORM( J )
511: GROW = GROW / XJ
512: 70 CONTINUE
513: END IF
514: 80 CONTINUE
515: END IF
516: *
517: IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
518: *
519: * Use the Level 2 BLAS solve if the reciprocal of the bound on
520: * elements of X is not too small.
521: *
522: CALL DTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
523: ELSE
524: *
525: * Use a Level 1 BLAS solve, scaling intermediate results.
526: *
527: IF( XMAX.GT.BIGNUM ) THEN
528: *
529: * Scale X so that its components are less than or equal to
530: * BIGNUM in absolute value.
531: *
532: SCALE = BIGNUM / XMAX
533: CALL DSCAL( N, SCALE, X, 1 )
534: XMAX = BIGNUM
535: END IF
536: *
537: IF( NOTRAN ) THEN
538: *
539: * Solve A * x = b
540: *
541: DO 110 J = JFIRST, JLAST, JINC
542: *
543: * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
544: *
545: XJ = ABS( X( J ) )
546: IF( NOUNIT ) THEN
547: TJJS = AB( MAIND, J )*TSCAL
548: ELSE
549: TJJS = TSCAL
550: IF( TSCAL.EQ.ONE )
551: $ GO TO 100
552: END IF
553: TJJ = ABS( TJJS )
554: IF( TJJ.GT.SMLNUM ) THEN
555: *
556: * abs(A(j,j)) > SMLNUM:
557: *
558: IF( TJJ.LT.ONE ) THEN
559: IF( XJ.GT.TJJ*BIGNUM ) THEN
560: *
561: * Scale x by 1/b(j).
562: *
563: REC = ONE / XJ
564: CALL DSCAL( N, REC, X, 1 )
565: SCALE = SCALE*REC
566: XMAX = XMAX*REC
567: END IF
568: END IF
569: X( J ) = X( J ) / TJJS
570: XJ = ABS( X( J ) )
571: ELSE IF( TJJ.GT.ZERO ) THEN
572: *
573: * 0 < abs(A(j,j)) <= SMLNUM:
574: *
575: IF( XJ.GT.TJJ*BIGNUM ) THEN
576: *
577: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
578: * to avoid overflow when dividing by A(j,j).
579: *
580: REC = ( TJJ*BIGNUM ) / XJ
581: IF( CNORM( J ).GT.ONE ) THEN
582: *
583: * Scale by 1/CNORM(j) to avoid overflow when
584: * multiplying x(j) times column j.
585: *
586: REC = REC / CNORM( J )
587: END IF
588: CALL DSCAL( N, REC, X, 1 )
589: SCALE = SCALE*REC
590: XMAX = XMAX*REC
591: END IF
592: X( J ) = X( J ) / TJJS
593: XJ = ABS( X( J ) )
594: ELSE
595: *
596: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
597: * scale = 0, and compute a solution to A*x = 0.
598: *
599: DO 90 I = 1, N
600: X( I ) = ZERO
601: 90 CONTINUE
602: X( J ) = ONE
603: XJ = ONE
604: SCALE = ZERO
605: XMAX = ZERO
606: END IF
607: 100 CONTINUE
608: *
609: * Scale x if necessary to avoid overflow when adding a
610: * multiple of column j of A.
611: *
612: IF( XJ.GT.ONE ) THEN
613: REC = ONE / XJ
614: IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
615: *
616: * Scale x by 1/(2*abs(x(j))).
617: *
618: REC = REC*HALF
619: CALL DSCAL( N, REC, X, 1 )
620: SCALE = SCALE*REC
621: END IF
622: ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
623: *
624: * Scale x by 1/2.
625: *
626: CALL DSCAL( N, HALF, X, 1 )
627: SCALE = SCALE*HALF
628: END IF
629: *
630: IF( UPPER ) THEN
631: IF( J.GT.1 ) THEN
632: *
633: * Compute the update
634: * x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
635: * x(j)* A(max(1,j-kd):j-1,j)
636: *
637: JLEN = MIN( KD, J-1 )
638: CALL DAXPY( JLEN, -X( J )*TSCAL,
639: $ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
640: I = IDAMAX( J-1, X, 1 )
641: XMAX = ABS( X( I ) )
642: END IF
643: ELSE IF( J.LT.N ) THEN
644: *
645: * Compute the update
646: * x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
647: * x(j) * A(j+1:min(j+kd,n),j)
648: *
649: JLEN = MIN( KD, N-J )
650: IF( JLEN.GT.0 )
651: $ CALL DAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
652: $ X( J+1 ), 1 )
653: I = J + IDAMAX( N-J, X( J+1 ), 1 )
654: XMAX = ABS( X( I ) )
655: END IF
656: 110 CONTINUE
657: *
658: ELSE
659: *
1.8 bertrand 660: * Solve A**T * x = b
1.1 bertrand 661: *
662: DO 160 J = JFIRST, JLAST, JINC
663: *
664: * Compute x(j) = b(j) - sum A(k,j)*x(k).
665: * k<>j
666: *
667: XJ = ABS( X( J ) )
668: USCAL = TSCAL
669: REC = ONE / MAX( XMAX, ONE )
670: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
671: *
672: * If x(j) could overflow, scale x by 1/(2*XMAX).
673: *
674: REC = REC*HALF
675: IF( NOUNIT ) THEN
676: TJJS = AB( MAIND, J )*TSCAL
677: ELSE
678: TJJS = TSCAL
679: END IF
680: TJJ = ABS( TJJS )
681: IF( TJJ.GT.ONE ) THEN
682: *
683: * Divide by A(j,j) when scaling x if A(j,j) > 1.
684: *
685: REC = MIN( ONE, REC*TJJ )
686: USCAL = USCAL / TJJS
687: END IF
688: IF( REC.LT.ONE ) THEN
689: CALL DSCAL( N, REC, X, 1 )
690: SCALE = SCALE*REC
691: XMAX = XMAX*REC
692: END IF
693: END IF
694: *
695: SUMJ = ZERO
696: IF( USCAL.EQ.ONE ) THEN
697: *
698: * If the scaling needed for A in the dot product is 1,
699: * call DDOT to perform the dot product.
700: *
701: IF( UPPER ) THEN
702: JLEN = MIN( KD, J-1 )
703: SUMJ = DDOT( JLEN, AB( KD+1-JLEN, J ), 1,
704: $ X( J-JLEN ), 1 )
705: ELSE
706: JLEN = MIN( KD, N-J )
707: IF( JLEN.GT.0 )
708: $ SUMJ = DDOT( JLEN, AB( 2, J ), 1, X( J+1 ), 1 )
709: END IF
710: ELSE
711: *
712: * Otherwise, use in-line code for the dot product.
713: *
714: IF( UPPER ) THEN
715: JLEN = MIN( KD, J-1 )
716: DO 120 I = 1, JLEN
717: SUMJ = SUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
718: $ X( J-JLEN-1+I )
719: 120 CONTINUE
720: ELSE
721: JLEN = MIN( KD, N-J )
722: DO 130 I = 1, JLEN
723: SUMJ = SUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
724: 130 CONTINUE
725: END IF
726: END IF
727: *
728: IF( USCAL.EQ.TSCAL ) THEN
729: *
730: * Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
731: * was not used to scale the dotproduct.
732: *
733: X( J ) = X( J ) - SUMJ
734: XJ = ABS( X( J ) )
735: IF( NOUNIT ) THEN
736: *
737: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
738: *
739: TJJS = AB( MAIND, J )*TSCAL
740: ELSE
741: TJJS = TSCAL
742: IF( TSCAL.EQ.ONE )
743: $ GO TO 150
744: END IF
745: TJJ = ABS( TJJS )
746: IF( TJJ.GT.SMLNUM ) THEN
747: *
748: * abs(A(j,j)) > SMLNUM:
749: *
750: IF( TJJ.LT.ONE ) THEN
751: IF( XJ.GT.TJJ*BIGNUM ) THEN
752: *
753: * Scale X by 1/abs(x(j)).
754: *
755: REC = ONE / XJ
756: CALL DSCAL( N, REC, X, 1 )
757: SCALE = SCALE*REC
758: XMAX = XMAX*REC
759: END IF
760: END IF
761: X( J ) = X( J ) / TJJS
762: ELSE IF( TJJ.GT.ZERO ) THEN
763: *
764: * 0 < abs(A(j,j)) <= SMLNUM:
765: *
766: IF( XJ.GT.TJJ*BIGNUM ) THEN
767: *
768: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
769: *
770: REC = ( TJJ*BIGNUM ) / XJ
771: CALL DSCAL( N, REC, X, 1 )
772: SCALE = SCALE*REC
773: XMAX = XMAX*REC
774: END IF
775: X( J ) = X( J ) / TJJS
776: ELSE
777: *
778: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
1.8 bertrand 779: * scale = 0, and compute a solution to A**T*x = 0.
1.1 bertrand 780: *
781: DO 140 I = 1, N
782: X( I ) = ZERO
783: 140 CONTINUE
784: X( J ) = ONE
785: SCALE = ZERO
786: XMAX = ZERO
787: END IF
788: 150 CONTINUE
789: ELSE
790: *
791: * Compute x(j) := x(j) / A(j,j) - sumj if the dot
792: * product has already been divided by 1/A(j,j).
793: *
794: X( J ) = X( J ) / TJJS - SUMJ
795: END IF
796: XMAX = MAX( XMAX, ABS( X( J ) ) )
797: 160 CONTINUE
798: END IF
799: SCALE = SCALE / TSCAL
800: END IF
801: *
802: * Scale the column norms by 1/TSCAL for return.
803: *
804: IF( TSCAL.NE.ONE ) THEN
805: CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
806: END IF
807: *
808: RETURN
809: *
810: * End of DLATBS
811: *
812: END