File:
[local] /
rpl /
lapack /
lapack /
zlatbs.f
Revision
1.15:
download - view:
text,
annotated -
select for diffs -
revision graph
Mon Jan 27 09:28:40 2014 UTC (10 years, 3 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_24,
rpl-4_1_23,
rpl-4_1_22,
rpl-4_1_21,
rpl-4_1_20,
rpl-4_1_19,
rpl-4_1_18,
rpl-4_1_17,
HEAD
Cohérence.
1: *> \brief \b ZLATBS 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 ZLATBS + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlatbs.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlatbs.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlatbs.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLATBS( 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 CNORM( * )
31: * COMPLEX*16 AB( LDAB, * ), X( * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> ZLATBS 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 ZTBSV 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*16 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*16 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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: *> \date September 2012
167: *
168: *> \ingroup complex16OTHERauxiliary
169: *
170: *> \par Further Details:
171: * =====================
172: *>
173: *> \verbatim
174: *>
175: *> A rough bound on x is computed; if that is less than overflow, ZTBSV
176: *> is called, otherwise, specific code is used which checks for possible
177: *> overflow or divide-by-zero at every operation.
178: *>
179: *> A columnwise scheme is used for solving A*x = b. The basic algorithm
180: *> if A is lower triangular is
181: *>
182: *> x[1:n] := b[1:n]
183: *> for j = 1, ..., n
184: *> x(j) := x(j) / A(j,j)
185: *> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
186: *> end
187: *>
188: *> Define bounds on the components of x after j iterations of the loop:
189: *> M(j) = bound on x[1:j]
190: *> G(j) = bound on x[j+1:n]
191: *> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
192: *>
193: *> Then for iteration j+1 we have
194: *> M(j+1) <= G(j) / | A(j+1,j+1) |
195: *> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
196: *> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
197: *>
198: *> where CNORM(j+1) is greater than or equal to the infinity-norm of
199: *> column j+1 of A, not counting the diagonal. Hence
200: *>
201: *> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
202: *> 1<=i<=j
203: *> and
204: *>
205: *> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
206: *> 1<=i< j
207: *>
208: *> Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTBSV if the
209: *> reciprocal of the largest M(j), j=1,..,n, is larger than
210: *> max(underflow, 1/overflow).
211: *>
212: *> The bound on x(j) is also used to determine when a step in the
213: *> columnwise method can be performed without fear of overflow. If
214: *> the computed bound is greater than a large constant, x is scaled to
215: *> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
216: *> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
217: *>
218: *> Similarly, a row-wise scheme is used to solve A**T *x = b or
219: *> A**H *x = b. The basic algorithm for A upper triangular is
220: *>
221: *> for j = 1, ..., n
222: *> x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
223: *> end
224: *>
225: *> We simultaneously compute two bounds
226: *> G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
227: *> M(j) = bound on x(i), 1<=i<=j
228: *>
229: *> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
230: *> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
231: *> Then the bound on x(j) is
232: *>
233: *> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
234: *>
235: *> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
236: *> 1<=i<=j
237: *>
238: *> and we can safely call ZTBSV if 1/M(n) and 1/G(n) are both greater
239: *> than max(underflow, 1/overflow).
240: *> \endverbatim
241: *>
242: * =====================================================================
243: SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
244: $ SCALE, CNORM, INFO )
245: *
246: * -- LAPACK auxiliary routine (version 3.4.2) --
247: * -- LAPACK is a software package provided by Univ. of Tennessee, --
248: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
249: * September 2012
250: *
251: * .. Scalar Arguments ..
252: CHARACTER DIAG, NORMIN, TRANS, UPLO
253: INTEGER INFO, KD, LDAB, N
254: DOUBLE PRECISION SCALE
255: * ..
256: * .. Array Arguments ..
257: DOUBLE PRECISION CNORM( * )
258: COMPLEX*16 AB( LDAB, * ), X( * )
259: * ..
260: *
261: * =====================================================================
262: *
263: * .. Parameters ..
264: DOUBLE PRECISION ZERO, HALF, ONE, TWO
265: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
266: $ TWO = 2.0D+0 )
267: * ..
268: * .. Local Scalars ..
269: LOGICAL NOTRAN, NOUNIT, UPPER
270: INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
271: DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
272: $ XBND, XJ, XMAX
273: COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
274: * ..
275: * .. External Functions ..
276: LOGICAL LSAME
277: INTEGER IDAMAX, IZAMAX
278: DOUBLE PRECISION DLAMCH, DZASUM
279: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
280: EXTERNAL LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC,
281: $ ZDOTU, ZLADIV
282: * ..
283: * .. External Subroutines ..
284: EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV
285: * ..
286: * .. Intrinsic Functions ..
287: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
288: * ..
289: * .. Statement Functions ..
290: DOUBLE PRECISION CABS1, CABS2
291: * ..
292: * .. Statement Function definitions ..
293: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
294: CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) +
295: $ ABS( DIMAG( ZDUM ) / 2.D0 )
296: * ..
297: * .. Executable Statements ..
298: *
299: INFO = 0
300: UPPER = LSAME( UPLO, 'U' )
301: NOTRAN = LSAME( TRANS, 'N' )
302: NOUNIT = LSAME( DIAG, 'N' )
303: *
304: * Test the input parameters.
305: *
306: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
307: INFO = -1
308: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
309: $ LSAME( TRANS, 'C' ) ) THEN
310: INFO = -2
311: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
312: INFO = -3
313: ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
314: $ LSAME( NORMIN, 'N' ) ) THEN
315: INFO = -4
316: ELSE IF( N.LT.0 ) THEN
317: INFO = -5
318: ELSE IF( KD.LT.0 ) THEN
319: INFO = -6
320: ELSE IF( LDAB.LT.KD+1 ) THEN
321: INFO = -8
322: END IF
323: IF( INFO.NE.0 ) THEN
324: CALL XERBLA( 'ZLATBS', -INFO )
325: RETURN
326: END IF
327: *
328: * Quick return if possible
329: *
330: IF( N.EQ.0 )
331: $ RETURN
332: *
333: * Determine machine dependent parameters to control overflow.
334: *
335: SMLNUM = DLAMCH( 'Safe minimum' )
336: BIGNUM = ONE / SMLNUM
337: CALL DLABAD( SMLNUM, BIGNUM )
338: SMLNUM = SMLNUM / DLAMCH( 'Precision' )
339: BIGNUM = ONE / SMLNUM
340: SCALE = ONE
341: *
342: IF( LSAME( NORMIN, 'N' ) ) THEN
343: *
344: * Compute the 1-norm of each column, not including the diagonal.
345: *
346: IF( UPPER ) THEN
347: *
348: * A is upper triangular.
349: *
350: DO 10 J = 1, N
351: JLEN = MIN( KD, J-1 )
352: CNORM( J ) = DZASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
353: 10 CONTINUE
354: ELSE
355: *
356: * A is lower triangular.
357: *
358: DO 20 J = 1, N
359: JLEN = MIN( KD, N-J )
360: IF( JLEN.GT.0 ) THEN
361: CNORM( J ) = DZASUM( JLEN, AB( 2, J ), 1 )
362: ELSE
363: CNORM( J ) = ZERO
364: END IF
365: 20 CONTINUE
366: END IF
367: END IF
368: *
369: * Scale the column norms by TSCAL if the maximum element in CNORM is
370: * greater than BIGNUM/2.
371: *
372: IMAX = IDAMAX( N, CNORM, 1 )
373: TMAX = CNORM( IMAX )
374: IF( TMAX.LE.BIGNUM*HALF ) THEN
375: TSCAL = ONE
376: ELSE
377: TSCAL = HALF / ( SMLNUM*TMAX )
378: CALL DSCAL( N, TSCAL, CNORM, 1 )
379: END IF
380: *
381: * Compute a bound on the computed solution vector to see if the
382: * Level 2 BLAS routine ZTBSV can be used.
383: *
384: XMAX = ZERO
385: DO 30 J = 1, N
386: XMAX = MAX( XMAX, CABS2( X( J ) ) )
387: 30 CONTINUE
388: XBND = XMAX
389: IF( NOTRAN ) THEN
390: *
391: * Compute the growth in A * x = b.
392: *
393: IF( UPPER ) THEN
394: JFIRST = N
395: JLAST = 1
396: JINC = -1
397: MAIND = KD + 1
398: ELSE
399: JFIRST = 1
400: JLAST = N
401: JINC = 1
402: MAIND = 1
403: END IF
404: *
405: IF( TSCAL.NE.ONE ) THEN
406: GROW = ZERO
407: GO TO 60
408: END IF
409: *
410: IF( NOUNIT ) THEN
411: *
412: * A is non-unit triangular.
413: *
414: * Compute GROW = 1/G(j) and XBND = 1/M(j).
415: * Initially, G(0) = max{x(i), i=1,...,n}.
416: *
417: GROW = HALF / MAX( XBND, SMLNUM )
418: XBND = GROW
419: DO 40 J = JFIRST, JLAST, JINC
420: *
421: * Exit the loop if the growth factor is too small.
422: *
423: IF( GROW.LE.SMLNUM )
424: $ GO TO 60
425: *
426: TJJS = AB( MAIND, J )
427: TJJ = CABS1( TJJS )
428: *
429: IF( TJJ.GE.SMLNUM ) THEN
430: *
431: * M(j) = G(j-1) / abs(A(j,j))
432: *
433: XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
434: ELSE
435: *
436: * M(j) could overflow, set XBND to 0.
437: *
438: XBND = ZERO
439: END IF
440: *
441: IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
442: *
443: * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
444: *
445: GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
446: ELSE
447: *
448: * G(j) could overflow, set GROW to 0.
449: *
450: GROW = ZERO
451: END IF
452: 40 CONTINUE
453: GROW = XBND
454: ELSE
455: *
456: * A is unit triangular.
457: *
458: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
459: *
460: GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
461: DO 50 J = JFIRST, JLAST, JINC
462: *
463: * Exit the loop if the growth factor is too small.
464: *
465: IF( GROW.LE.SMLNUM )
466: $ GO TO 60
467: *
468: * G(j) = G(j-1)*( 1 + CNORM(j) )
469: *
470: GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
471: 50 CONTINUE
472: END IF
473: 60 CONTINUE
474: *
475: ELSE
476: *
477: * Compute the growth in A**T * x = b or A**H * x = b.
478: *
479: IF( UPPER ) THEN
480: JFIRST = 1
481: JLAST = N
482: JINC = 1
483: MAIND = KD + 1
484: ELSE
485: JFIRST = N
486: JLAST = 1
487: JINC = -1
488: MAIND = 1
489: END IF
490: *
491: IF( TSCAL.NE.ONE ) THEN
492: GROW = ZERO
493: GO TO 90
494: END IF
495: *
496: IF( NOUNIT ) THEN
497: *
498: * A is non-unit triangular.
499: *
500: * Compute GROW = 1/G(j) and XBND = 1/M(j).
501: * Initially, M(0) = max{x(i), i=1,...,n}.
502: *
503: GROW = HALF / MAX( XBND, SMLNUM )
504: XBND = GROW
505: DO 70 J = JFIRST, JLAST, JINC
506: *
507: * Exit the loop if the growth factor is too small.
508: *
509: IF( GROW.LE.SMLNUM )
510: $ GO TO 90
511: *
512: * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
513: *
514: XJ = ONE + CNORM( J )
515: GROW = MIN( GROW, XBND / XJ )
516: *
517: TJJS = AB( MAIND, J )
518: TJJ = CABS1( TJJS )
519: *
520: IF( TJJ.GE.SMLNUM ) THEN
521: *
522: * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
523: *
524: IF( XJ.GT.TJJ )
525: $ XBND = XBND*( TJJ / XJ )
526: ELSE
527: *
528: * M(j) could overflow, set XBND to 0.
529: *
530: XBND = ZERO
531: END IF
532: 70 CONTINUE
533: GROW = MIN( GROW, XBND )
534: ELSE
535: *
536: * A is unit triangular.
537: *
538: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
539: *
540: GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
541: DO 80 J = JFIRST, JLAST, JINC
542: *
543: * Exit the loop if the growth factor is too small.
544: *
545: IF( GROW.LE.SMLNUM )
546: $ GO TO 90
547: *
548: * G(j) = ( 1 + CNORM(j) )*G(j-1)
549: *
550: XJ = ONE + CNORM( J )
551: GROW = GROW / XJ
552: 80 CONTINUE
553: END IF
554: 90 CONTINUE
555: END IF
556: *
557: IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
558: *
559: * Use the Level 2 BLAS solve if the reciprocal of the bound on
560: * elements of X is not too small.
561: *
562: CALL ZTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
563: ELSE
564: *
565: * Use a Level 1 BLAS solve, scaling intermediate results.
566: *
567: IF( XMAX.GT.BIGNUM*HALF ) THEN
568: *
569: * Scale X so that its components are less than or equal to
570: * BIGNUM in absolute value.
571: *
572: SCALE = ( BIGNUM*HALF ) / XMAX
573: CALL ZDSCAL( N, SCALE, X, 1 )
574: XMAX = BIGNUM
575: ELSE
576: XMAX = XMAX*TWO
577: END IF
578: *
579: IF( NOTRAN ) THEN
580: *
581: * Solve A * x = b
582: *
583: DO 120 J = JFIRST, JLAST, JINC
584: *
585: * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
586: *
587: XJ = CABS1( X( J ) )
588: IF( NOUNIT ) THEN
589: TJJS = AB( MAIND, J )*TSCAL
590: ELSE
591: TJJS = TSCAL
592: IF( TSCAL.EQ.ONE )
593: $ GO TO 110
594: END IF
595: TJJ = CABS1( TJJS )
596: IF( TJJ.GT.SMLNUM ) THEN
597: *
598: * abs(A(j,j)) > SMLNUM:
599: *
600: IF( TJJ.LT.ONE ) THEN
601: IF( XJ.GT.TJJ*BIGNUM ) THEN
602: *
603: * Scale x by 1/b(j).
604: *
605: REC = ONE / XJ
606: CALL ZDSCAL( N, REC, X, 1 )
607: SCALE = SCALE*REC
608: XMAX = XMAX*REC
609: END IF
610: END IF
611: X( J ) = ZLADIV( X( J ), TJJS )
612: XJ = CABS1( X( J ) )
613: ELSE IF( TJJ.GT.ZERO ) THEN
614: *
615: * 0 < abs(A(j,j)) <= SMLNUM:
616: *
617: IF( XJ.GT.TJJ*BIGNUM ) THEN
618: *
619: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
620: * to avoid overflow when dividing by A(j,j).
621: *
622: REC = ( TJJ*BIGNUM ) / XJ
623: IF( CNORM( J ).GT.ONE ) THEN
624: *
625: * Scale by 1/CNORM(j) to avoid overflow when
626: * multiplying x(j) times column j.
627: *
628: REC = REC / CNORM( J )
629: END IF
630: CALL ZDSCAL( N, REC, X, 1 )
631: SCALE = SCALE*REC
632: XMAX = XMAX*REC
633: END IF
634: X( J ) = ZLADIV( X( J ), TJJS )
635: XJ = CABS1( X( J ) )
636: ELSE
637: *
638: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
639: * scale = 0, and compute a solution to A*x = 0.
640: *
641: DO 100 I = 1, N
642: X( I ) = ZERO
643: 100 CONTINUE
644: X( J ) = ONE
645: XJ = ONE
646: SCALE = ZERO
647: XMAX = ZERO
648: END IF
649: 110 CONTINUE
650: *
651: * Scale x if necessary to avoid overflow when adding a
652: * multiple of column j of A.
653: *
654: IF( XJ.GT.ONE ) THEN
655: REC = ONE / XJ
656: IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
657: *
658: * Scale x by 1/(2*abs(x(j))).
659: *
660: REC = REC*HALF
661: CALL ZDSCAL( N, REC, X, 1 )
662: SCALE = SCALE*REC
663: END IF
664: ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
665: *
666: * Scale x by 1/2.
667: *
668: CALL ZDSCAL( N, HALF, X, 1 )
669: SCALE = SCALE*HALF
670: END IF
671: *
672: IF( UPPER ) THEN
673: IF( J.GT.1 ) THEN
674: *
675: * Compute the update
676: * x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
677: * x(j)* A(max(1,j-kd):j-1,j)
678: *
679: JLEN = MIN( KD, J-1 )
680: CALL ZAXPY( JLEN, -X( J )*TSCAL,
681: $ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
682: I = IZAMAX( J-1, X, 1 )
683: XMAX = CABS1( X( I ) )
684: END IF
685: ELSE IF( J.LT.N ) THEN
686: *
687: * Compute the update
688: * x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
689: * x(j) * A(j+1:min(j+kd,n),j)
690: *
691: JLEN = MIN( KD, N-J )
692: IF( JLEN.GT.0 )
693: $ CALL ZAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
694: $ X( J+1 ), 1 )
695: I = J + IZAMAX( N-J, X( J+1 ), 1 )
696: XMAX = CABS1( X( I ) )
697: END IF
698: 120 CONTINUE
699: *
700: ELSE IF( LSAME( TRANS, 'T' ) ) THEN
701: *
702: * Solve A**T * x = b
703: *
704: DO 170 J = JFIRST, JLAST, JINC
705: *
706: * Compute x(j) = b(j) - sum A(k,j)*x(k).
707: * k<>j
708: *
709: XJ = CABS1( X( J ) )
710: USCAL = TSCAL
711: REC = ONE / MAX( XMAX, ONE )
712: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
713: *
714: * If x(j) could overflow, scale x by 1/(2*XMAX).
715: *
716: REC = REC*HALF
717: IF( NOUNIT ) THEN
718: TJJS = AB( MAIND, J )*TSCAL
719: ELSE
720: TJJS = TSCAL
721: END IF
722: TJJ = CABS1( TJJS )
723: IF( TJJ.GT.ONE ) THEN
724: *
725: * Divide by A(j,j) when scaling x if A(j,j) > 1.
726: *
727: REC = MIN( ONE, REC*TJJ )
728: USCAL = ZLADIV( USCAL, TJJS )
729: END IF
730: IF( REC.LT.ONE ) THEN
731: CALL ZDSCAL( N, REC, X, 1 )
732: SCALE = SCALE*REC
733: XMAX = XMAX*REC
734: END IF
735: END IF
736: *
737: CSUMJ = ZERO
738: IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
739: *
740: * If the scaling needed for A in the dot product is 1,
741: * call ZDOTU to perform the dot product.
742: *
743: IF( UPPER ) THEN
744: JLEN = MIN( KD, J-1 )
745: CSUMJ = ZDOTU( JLEN, AB( KD+1-JLEN, J ), 1,
746: $ X( J-JLEN ), 1 )
747: ELSE
748: JLEN = MIN( KD, N-J )
749: IF( JLEN.GT.1 )
750: $ CSUMJ = ZDOTU( JLEN, AB( 2, J ), 1, X( J+1 ),
751: $ 1 )
752: END IF
753: ELSE
754: *
755: * Otherwise, use in-line code for the dot product.
756: *
757: IF( UPPER ) THEN
758: JLEN = MIN( KD, J-1 )
759: DO 130 I = 1, JLEN
760: CSUMJ = CSUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
761: $ X( J-JLEN-1+I )
762: 130 CONTINUE
763: ELSE
764: JLEN = MIN( KD, N-J )
765: DO 140 I = 1, JLEN
766: CSUMJ = CSUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
767: 140 CONTINUE
768: END IF
769: END IF
770: *
771: IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
772: *
773: * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
774: * was not used to scale the dotproduct.
775: *
776: X( J ) = X( J ) - CSUMJ
777: XJ = CABS1( X( J ) )
778: IF( NOUNIT ) THEN
779: *
780: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
781: *
782: TJJS = AB( MAIND, J )*TSCAL
783: ELSE
784: TJJS = TSCAL
785: IF( TSCAL.EQ.ONE )
786: $ GO TO 160
787: END IF
788: TJJ = CABS1( TJJS )
789: IF( TJJ.GT.SMLNUM ) THEN
790: *
791: * abs(A(j,j)) > SMLNUM:
792: *
793: IF( TJJ.LT.ONE ) THEN
794: IF( XJ.GT.TJJ*BIGNUM ) THEN
795: *
796: * Scale X by 1/abs(x(j)).
797: *
798: REC = ONE / XJ
799: CALL ZDSCAL( N, REC, X, 1 )
800: SCALE = SCALE*REC
801: XMAX = XMAX*REC
802: END IF
803: END IF
804: X( J ) = ZLADIV( X( J ), TJJS )
805: ELSE IF( TJJ.GT.ZERO ) THEN
806: *
807: * 0 < abs(A(j,j)) <= SMLNUM:
808: *
809: IF( XJ.GT.TJJ*BIGNUM ) THEN
810: *
811: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
812: *
813: REC = ( TJJ*BIGNUM ) / XJ
814: CALL ZDSCAL( N, REC, X, 1 )
815: SCALE = SCALE*REC
816: XMAX = XMAX*REC
817: END IF
818: X( J ) = ZLADIV( X( J ), TJJS )
819: ELSE
820: *
821: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
822: * scale = 0 and compute a solution to A**T *x = 0.
823: *
824: DO 150 I = 1, N
825: X( I ) = ZERO
826: 150 CONTINUE
827: X( J ) = ONE
828: SCALE = ZERO
829: XMAX = ZERO
830: END IF
831: 160 CONTINUE
832: ELSE
833: *
834: * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
835: * product has already been divided by 1/A(j,j).
836: *
837: X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
838: END IF
839: XMAX = MAX( XMAX, CABS1( X( J ) ) )
840: 170 CONTINUE
841: *
842: ELSE
843: *
844: * Solve A**H * x = b
845: *
846: DO 220 J = JFIRST, JLAST, JINC
847: *
848: * Compute x(j) = b(j) - sum A(k,j)*x(k).
849: * k<>j
850: *
851: XJ = CABS1( X( J ) )
852: USCAL = TSCAL
853: REC = ONE / MAX( XMAX, ONE )
854: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
855: *
856: * If x(j) could overflow, scale x by 1/(2*XMAX).
857: *
858: REC = REC*HALF
859: IF( NOUNIT ) THEN
860: TJJS = DCONJG( AB( MAIND, J ) )*TSCAL
861: ELSE
862: TJJS = TSCAL
863: END IF
864: TJJ = CABS1( TJJS )
865: IF( TJJ.GT.ONE ) THEN
866: *
867: * Divide by A(j,j) when scaling x if A(j,j) > 1.
868: *
869: REC = MIN( ONE, REC*TJJ )
870: USCAL = ZLADIV( USCAL, TJJS )
871: END IF
872: IF( REC.LT.ONE ) THEN
873: CALL ZDSCAL( N, REC, X, 1 )
874: SCALE = SCALE*REC
875: XMAX = XMAX*REC
876: END IF
877: END IF
878: *
879: CSUMJ = ZERO
880: IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
881: *
882: * If the scaling needed for A in the dot product is 1,
883: * call ZDOTC to perform the dot product.
884: *
885: IF( UPPER ) THEN
886: JLEN = MIN( KD, J-1 )
887: CSUMJ = ZDOTC( JLEN, AB( KD+1-JLEN, J ), 1,
888: $ X( J-JLEN ), 1 )
889: ELSE
890: JLEN = MIN( KD, N-J )
891: IF( JLEN.GT.1 )
892: $ CSUMJ = ZDOTC( JLEN, AB( 2, J ), 1, X( J+1 ),
893: $ 1 )
894: END IF
895: ELSE
896: *
897: * Otherwise, use in-line code for the dot product.
898: *
899: IF( UPPER ) THEN
900: JLEN = MIN( KD, J-1 )
901: DO 180 I = 1, JLEN
902: CSUMJ = CSUMJ + ( DCONJG( AB( KD+I-JLEN, J ) )*
903: $ USCAL )*X( J-JLEN-1+I )
904: 180 CONTINUE
905: ELSE
906: JLEN = MIN( KD, N-J )
907: DO 190 I = 1, JLEN
908: CSUMJ = CSUMJ + ( DCONJG( AB( I+1, J ) )*USCAL )
909: $ *X( J+I )
910: 190 CONTINUE
911: END IF
912: END IF
913: *
914: IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
915: *
916: * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
917: * was not used to scale the dotproduct.
918: *
919: X( J ) = X( J ) - CSUMJ
920: XJ = CABS1( X( J ) )
921: IF( NOUNIT ) THEN
922: *
923: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
924: *
925: TJJS = DCONJG( AB( MAIND, J ) )*TSCAL
926: ELSE
927: TJJS = TSCAL
928: IF( TSCAL.EQ.ONE )
929: $ GO TO 210
930: END IF
931: TJJ = CABS1( TJJS )
932: IF( TJJ.GT.SMLNUM ) THEN
933: *
934: * abs(A(j,j)) > SMLNUM:
935: *
936: IF( TJJ.LT.ONE ) THEN
937: IF( XJ.GT.TJJ*BIGNUM ) THEN
938: *
939: * Scale X by 1/abs(x(j)).
940: *
941: REC = ONE / XJ
942: CALL ZDSCAL( N, REC, X, 1 )
943: SCALE = SCALE*REC
944: XMAX = XMAX*REC
945: END IF
946: END IF
947: X( J ) = ZLADIV( X( J ), TJJS )
948: ELSE IF( TJJ.GT.ZERO ) THEN
949: *
950: * 0 < abs(A(j,j)) <= SMLNUM:
951: *
952: IF( XJ.GT.TJJ*BIGNUM ) THEN
953: *
954: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
955: *
956: REC = ( TJJ*BIGNUM ) / XJ
957: CALL ZDSCAL( N, REC, X, 1 )
958: SCALE = SCALE*REC
959: XMAX = XMAX*REC
960: END IF
961: X( J ) = ZLADIV( X( J ), TJJS )
962: ELSE
963: *
964: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
965: * scale = 0 and compute a solution to A**H *x = 0.
966: *
967: DO 200 I = 1, N
968: X( I ) = ZERO
969: 200 CONTINUE
970: X( J ) = ONE
971: SCALE = ZERO
972: XMAX = ZERO
973: END IF
974: 210 CONTINUE
975: ELSE
976: *
977: * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
978: * product has already been divided by 1/A(j,j).
979: *
980: X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
981: END IF
982: XMAX = MAX( XMAX, CABS1( X( J ) ) )
983: 220 CONTINUE
984: END IF
985: SCALE = SCALE / TSCAL
986: END IF
987: *
988: * Scale the column norms by 1/TSCAL for return.
989: *
990: IF( TSCAL.NE.ONE ) THEN
991: CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
992: END IF
993: *
994: RETURN
995: *
996: * End of ZLATBS
997: *
998: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlatbs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack