1: *> \brief \b DLATBS 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 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
137: *> CNORM is DOUBLE PRECISION 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 doubleOTHERauxiliary
166: *
167: *> \par Further Details:
168: * =====================
169: *>
170: *> \verbatim
171: *>
172: *> A rough bound on x is computed; if that is less than overflow, DTBSV
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 DTBSV 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 DTBSV if 1/M(n) and 1/G(n) are both greater
236: *> than max(underflow, 1/overflow).
237: *> \endverbatim
238: *>
239: * =====================================================================
240: SUBROUTINE DLATBS( 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: DOUBLE PRECISION SCALE
251: * ..
252: * .. Array Arguments ..
253: DOUBLE PRECISION AB( LDAB, * ), CNORM( * ), X( * )
254: * ..
255: *
256: * =====================================================================
257: *
258: * .. Parameters ..
259: DOUBLE PRECISION ZERO, HALF, ONE
260: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
261: * ..
262: * .. Local Scalars ..
263: LOGICAL NOTRAN, NOUNIT, UPPER
264: INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
265: DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
266: $ TMAX, TSCAL, USCAL, XBND, XJ, XMAX
267: * ..
268: * .. External Functions ..
269: LOGICAL LSAME
270: INTEGER IDAMAX
271: DOUBLE PRECISION DASUM, DDOT, DLAMCH
272: EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
273: * ..
274: * .. External Subroutines ..
275: EXTERNAL DAXPY, DSCAL, DTBSV, 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( 'DLATBS', -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 = DLAMCH( 'Safe minimum' ) / DLAMCH( '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 ) = DASUM( 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 ) = DASUM( 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 = IDAMAX( 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 DSCAL( 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 DTBSV can be used.
363: *
364: J = IDAMAX( 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 DTBSV( 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 DSCAL( N, SCALE, X, 1 )
531: XMAX = BIGNUM
532: END IF
533: *
534: IF( NOTRAN ) THEN
535: *
536: * Solve A * x = b
537: *
538: DO 110 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 100
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 DSCAL( 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 DSCAL( 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: 100 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 DSCAL( 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 DSCAL( 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 DAXPY( JLEN, -X( J )*TSCAL,
636: $ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
637: I = IDAMAX( 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 DAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
649: $ X( J+1 ), 1 )
650: I = J + IDAMAX( N-J, X( J+1 ), 1 )
651: XMAX = ABS( X( I ) )
652: END IF
653: 110 CONTINUE
654: *
655: ELSE
656: *
657: * Solve A**T * x = b
658: *
659: DO 160 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 DSCAL( 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 DDOT to perform the dot product.
697: *
698: IF( UPPER ) THEN
699: JLEN = MIN( KD, J-1 )
700: SUMJ = DDOT( 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 = DDOT( 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 120 I = 1, JLEN
714: SUMJ = SUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
715: $ X( J-JLEN-1+I )
716: 120 CONTINUE
717: ELSE
718: JLEN = MIN( KD, N-J )
719: DO 130 I = 1, JLEN
720: SUMJ = SUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
721: 130 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 150
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 DSCAL( 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 DSCAL( 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 140 I = 1, N
779: X( I ) = ZERO
780: 140 CONTINUE
781: X( J ) = ONE
782: SCALE = ZERO
783: XMAX = ZERO
784: END IF
785: 150 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: 160 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 DSCAL( N, ONE / TSCAL, CNORM, 1 )
803: END IF
804: *
805: RETURN
806: *
807: * End of DLATBS
808: *
809: END
CVSweb interface <joel.bertrand@systella.fr>
![[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