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: *> \ingroup complex16OTHERauxiliary
167: *
168: *> \par Further Details:
169: * =====================
170: *>
171: *> \verbatim
172: *>
173: *> A rough bound on x is computed; if that is less than overflow, ZTBSV
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 ZTBSV 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 ZTBSV if 1/M(n) and 1/G(n) are both greater
237: *> than max(underflow, 1/overflow).
238: *> \endverbatim
239: *>
240: * =====================================================================
241: SUBROUTINE ZLATBS( 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: DOUBLE PRECISION SCALE
252: * ..
253: * .. Array Arguments ..
254: DOUBLE PRECISION CNORM( * )
255: COMPLEX*16 AB( LDAB, * ), X( * )
256: * ..
257: *
258: * =====================================================================
259: *
260: * .. Parameters ..
261: DOUBLE PRECISION ZERO, HALF, ONE, TWO
262: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
263: $ TWO = 2.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, TJJ, TMAX, TSCAL,
269: $ XBND, XJ, XMAX
270: COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
271: * ..
272: * .. External Functions ..
273: LOGICAL LSAME
274: INTEGER IDAMAX, IZAMAX
275: DOUBLE PRECISION DLAMCH, DZASUM
276: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
277: EXTERNAL LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC,
278: $ ZDOTU, ZLADIV
279: * ..
280: * .. External Subroutines ..
281: EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV
282: * ..
283: * .. Intrinsic Functions ..
284: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
285: * ..
286: * .. Statement Functions ..
287: DOUBLE PRECISION CABS1, CABS2
288: * ..
289: * .. Statement Function definitions ..
290: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
291: CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) +
292: $ ABS( DIMAG( ZDUM ) / 2.D0 )
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( 'ZLATBS', -INFO )
322: RETURN
323: END IF
324: *
325: * Quick return if possible
326: *
327: SCALE = ONE
328: IF( N.EQ.0 )
329: $ RETURN
330: *
331: * Determine machine dependent parameters to control overflow.
332: *
333: SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
334: BIGNUM = ONE / SMLNUM
335: *
336: IF( LSAME( NORMIN, 'N' ) ) THEN
337: *
338: * Compute the 1-norm of each column, not including the diagonal.
339: *
340: IF( UPPER ) THEN
341: *
342: * A is upper triangular.
343: *
344: DO 10 J = 1, N
345: JLEN = MIN( KD, J-1 )
346: CNORM( J ) = DZASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
347: 10 CONTINUE
348: ELSE
349: *
350: * A is lower triangular.
351: *
352: DO 20 J = 1, N
353: JLEN = MIN( KD, N-J )
354: IF( JLEN.GT.0 ) THEN
355: CNORM( J ) = DZASUM( JLEN, AB( 2, J ), 1 )
356: ELSE
357: CNORM( J ) = ZERO
358: END IF
359: 20 CONTINUE
360: END IF
361: END IF
362: *
363: * Scale the column norms by TSCAL if the maximum element in CNORM is
364: * greater than BIGNUM/2.
365: *
366: IMAX = IDAMAX( N, CNORM, 1 )
367: TMAX = CNORM( IMAX )
368: IF( TMAX.LE.BIGNUM*HALF ) THEN
369: TSCAL = ONE
370: ELSE
371: TSCAL = HALF / ( SMLNUM*TMAX )
372: CALL DSCAL( N, TSCAL, CNORM, 1 )
373: END IF
374: *
375: * Compute a bound on the computed solution vector to see if the
376: * Level 2 BLAS routine ZTBSV can be used.
377: *
378: XMAX = ZERO
379: DO 30 J = 1, N
380: XMAX = MAX( XMAX, CABS2( X( J ) ) )
381: 30 CONTINUE
382: XBND = XMAX
383: IF( NOTRAN ) THEN
384: *
385: * Compute the growth in A * x = b.
386: *
387: IF( UPPER ) THEN
388: JFIRST = N
389: JLAST = 1
390: JINC = -1
391: MAIND = KD + 1
392: ELSE
393: JFIRST = 1
394: JLAST = N
395: JINC = 1
396: MAIND = 1
397: END IF
398: *
399: IF( TSCAL.NE.ONE ) THEN
400: GROW = ZERO
401: GO TO 60
402: END IF
403: *
404: IF( NOUNIT ) THEN
405: *
406: * A is non-unit triangular.
407: *
408: * Compute GROW = 1/G(j) and XBND = 1/M(j).
409: * Initially, G(0) = max{x(i), i=1,...,n}.
410: *
411: GROW = HALF / MAX( XBND, SMLNUM )
412: XBND = GROW
413: DO 40 J = JFIRST, JLAST, JINC
414: *
415: * Exit the loop if the growth factor is too small.
416: *
417: IF( GROW.LE.SMLNUM )
418: $ GO TO 60
419: *
420: TJJS = AB( MAIND, J )
421: TJJ = CABS1( TJJS )
422: *
423: IF( TJJ.GE.SMLNUM ) THEN
424: *
425: * M(j) = G(j-1) / abs(A(j,j))
426: *
427: XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
428: ELSE
429: *
430: * M(j) could overflow, set XBND to 0.
431: *
432: XBND = ZERO
433: END IF
434: *
435: IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
436: *
437: * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
438: *
439: GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
440: ELSE
441: *
442: * G(j) could overflow, set GROW to 0.
443: *
444: GROW = ZERO
445: END IF
446: 40 CONTINUE
447: GROW = XBND
448: ELSE
449: *
450: * A is unit triangular.
451: *
452: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
453: *
454: GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
455: DO 50 J = JFIRST, JLAST, JINC
456: *
457: * Exit the loop if the growth factor is too small.
458: *
459: IF( GROW.LE.SMLNUM )
460: $ GO TO 60
461: *
462: * G(j) = G(j-1)*( 1 + CNORM(j) )
463: *
464: GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
465: 50 CONTINUE
466: END IF
467: 60 CONTINUE
468: *
469: ELSE
470: *
471: * Compute the growth in A**T * x = b or A**H * x = b.
472: *
473: IF( UPPER ) THEN
474: JFIRST = 1
475: JLAST = N
476: JINC = 1
477: MAIND = KD + 1
478: ELSE
479: JFIRST = N
480: JLAST = 1
481: JINC = -1
482: MAIND = 1
483: END IF
484: *
485: IF( TSCAL.NE.ONE ) THEN
486: GROW = ZERO
487: GO TO 90
488: END IF
489: *
490: IF( NOUNIT ) THEN
491: *
492: * A is non-unit triangular.
493: *
494: * Compute GROW = 1/G(j) and XBND = 1/M(j).
495: * Initially, M(0) = max{x(i), i=1,...,n}.
496: *
497: GROW = HALF / MAX( XBND, SMLNUM )
498: XBND = GROW
499: DO 70 J = JFIRST, JLAST, JINC
500: *
501: * Exit the loop if the growth factor is too small.
502: *
503: IF( GROW.LE.SMLNUM )
504: $ GO TO 90
505: *
506: * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
507: *
508: XJ = ONE + CNORM( J )
509: GROW = MIN( GROW, XBND / XJ )
510: *
511: TJJS = AB( MAIND, J )
512: TJJ = CABS1( TJJS )
513: *
514: IF( TJJ.GE.SMLNUM ) THEN
515: *
516: * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
517: *
518: IF( XJ.GT.TJJ )
519: $ XBND = XBND*( TJJ / XJ )
520: ELSE
521: *
522: * M(j) could overflow, set XBND to 0.
523: *
524: XBND = ZERO
525: END IF
526: 70 CONTINUE
527: GROW = MIN( GROW, XBND )
528: ELSE
529: *
530: * A is unit triangular.
531: *
532: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
533: *
534: GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
535: DO 80 J = JFIRST, JLAST, JINC
536: *
537: * Exit the loop if the growth factor is too small.
538: *
539: IF( GROW.LE.SMLNUM )
540: $ GO TO 90
541: *
542: * G(j) = ( 1 + CNORM(j) )*G(j-1)
543: *
544: XJ = ONE + CNORM( J )
545: GROW = GROW / XJ
546: 80 CONTINUE
547: END IF
548: 90 CONTINUE
549: END IF
550: *
551: IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
552: *
553: * Use the Level 2 BLAS solve if the reciprocal of the bound on
554: * elements of X is not too small.
555: *
556: CALL ZTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
557: ELSE
558: *
559: * Use a Level 1 BLAS solve, scaling intermediate results.
560: *
561: IF( XMAX.GT.BIGNUM*HALF ) THEN
562: *
563: * Scale X so that its components are less than or equal to
564: * BIGNUM in absolute value.
565: *
566: SCALE = ( BIGNUM*HALF ) / XMAX
567: CALL ZDSCAL( N, SCALE, X, 1 )
568: XMAX = BIGNUM
569: ELSE
570: XMAX = XMAX*TWO
571: END IF
572: *
573: IF( NOTRAN ) THEN
574: *
575: * Solve A * x = b
576: *
577: DO 120 J = JFIRST, JLAST, JINC
578: *
579: * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
580: *
581: XJ = CABS1( X( J ) )
582: IF( NOUNIT ) THEN
583: TJJS = AB( MAIND, J )*TSCAL
584: ELSE
585: TJJS = TSCAL
586: IF( TSCAL.EQ.ONE )
587: $ GO TO 110
588: END IF
589: TJJ = CABS1( TJJS )
590: IF( TJJ.GT.SMLNUM ) THEN
591: *
592: * abs(A(j,j)) > SMLNUM:
593: *
594: IF( TJJ.LT.ONE ) THEN
595: IF( XJ.GT.TJJ*BIGNUM ) THEN
596: *
597: * Scale x by 1/b(j).
598: *
599: REC = ONE / XJ
600: CALL ZDSCAL( N, REC, X, 1 )
601: SCALE = SCALE*REC
602: XMAX = XMAX*REC
603: END IF
604: END IF
605: X( J ) = ZLADIV( X( J ), TJJS )
606: XJ = CABS1( X( J ) )
607: ELSE IF( TJJ.GT.ZERO ) THEN
608: *
609: * 0 < abs(A(j,j)) <= SMLNUM:
610: *
611: IF( XJ.GT.TJJ*BIGNUM ) THEN
612: *
613: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
614: * to avoid overflow when dividing by A(j,j).
615: *
616: REC = ( TJJ*BIGNUM ) / XJ
617: IF( CNORM( J ).GT.ONE ) THEN
618: *
619: * Scale by 1/CNORM(j) to avoid overflow when
620: * multiplying x(j) times column j.
621: *
622: REC = REC / CNORM( J )
623: END IF
624: CALL ZDSCAL( N, REC, X, 1 )
625: SCALE = SCALE*REC
626: XMAX = XMAX*REC
627: END IF
628: X( J ) = ZLADIV( X( J ), TJJS )
629: XJ = CABS1( X( J ) )
630: ELSE
631: *
632: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
633: * scale = 0, and compute a solution to A*x = 0.
634: *
635: DO 100 I = 1, N
636: X( I ) = ZERO
637: 100 CONTINUE
638: X( J ) = ONE
639: XJ = ONE
640: SCALE = ZERO
641: XMAX = ZERO
642: END IF
643: 110 CONTINUE
644: *
645: * Scale x if necessary to avoid overflow when adding a
646: * multiple of column j of A.
647: *
648: IF( XJ.GT.ONE ) THEN
649: REC = ONE / XJ
650: IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
651: *
652: * Scale x by 1/(2*abs(x(j))).
653: *
654: REC = REC*HALF
655: CALL ZDSCAL( N, REC, X, 1 )
656: SCALE = SCALE*REC
657: END IF
658: ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
659: *
660: * Scale x by 1/2.
661: *
662: CALL ZDSCAL( N, HALF, X, 1 )
663: SCALE = SCALE*HALF
664: END IF
665: *
666: IF( UPPER ) THEN
667: IF( J.GT.1 ) THEN
668: *
669: * Compute the update
670: * x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
671: * x(j)* A(max(1,j-kd):j-1,j)
672: *
673: JLEN = MIN( KD, J-1 )
674: CALL ZAXPY( JLEN, -X( J )*TSCAL,
675: $ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
676: I = IZAMAX( J-1, X, 1 )
677: XMAX = CABS1( X( I ) )
678: END IF
679: ELSE IF( J.LT.N ) THEN
680: *
681: * Compute the update
682: * x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
683: * x(j) * A(j+1:min(j+kd,n),j)
684: *
685: JLEN = MIN( KD, N-J )
686: IF( JLEN.GT.0 )
687: $ CALL ZAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
688: $ X( J+1 ), 1 )
689: I = J + IZAMAX( N-J, X( J+1 ), 1 )
690: XMAX = CABS1( X( I ) )
691: END IF
692: 120 CONTINUE
693: *
694: ELSE IF( LSAME( TRANS, 'T' ) ) THEN
695: *
696: * Solve A**T * x = b
697: *
698: DO 170 J = JFIRST, JLAST, JINC
699: *
700: * Compute x(j) = b(j) - sum A(k,j)*x(k).
701: * k<>j
702: *
703: XJ = CABS1( X( J ) )
704: USCAL = TSCAL
705: REC = ONE / MAX( XMAX, ONE )
706: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
707: *
708: * If x(j) could overflow, scale x by 1/(2*XMAX).
709: *
710: REC = REC*HALF
711: IF( NOUNIT ) THEN
712: TJJS = AB( MAIND, J )*TSCAL
713: ELSE
714: TJJS = TSCAL
715: END IF
716: TJJ = CABS1( TJJS )
717: IF( TJJ.GT.ONE ) THEN
718: *
719: * Divide by A(j,j) when scaling x if A(j,j) > 1.
720: *
721: REC = MIN( ONE, REC*TJJ )
722: USCAL = ZLADIV( USCAL, TJJS )
723: END IF
724: IF( REC.LT.ONE ) THEN
725: CALL ZDSCAL( N, REC, X, 1 )
726: SCALE = SCALE*REC
727: XMAX = XMAX*REC
728: END IF
729: END IF
730: *
731: CSUMJ = ZERO
732: IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
733: *
734: * If the scaling needed for A in the dot product is 1,
735: * call ZDOTU to perform the dot product.
736: *
737: IF( UPPER ) THEN
738: JLEN = MIN( KD, J-1 )
739: CSUMJ = ZDOTU( JLEN, AB( KD+1-JLEN, J ), 1,
740: $ X( J-JLEN ), 1 )
741: ELSE
742: JLEN = MIN( KD, N-J )
743: IF( JLEN.GT.1 )
744: $ CSUMJ = ZDOTU( JLEN, AB( 2, J ), 1, X( J+1 ),
745: $ 1 )
746: END IF
747: ELSE
748: *
749: * Otherwise, use in-line code for the dot product.
750: *
751: IF( UPPER ) THEN
752: JLEN = MIN( KD, J-1 )
753: DO 130 I = 1, JLEN
754: CSUMJ = CSUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
755: $ X( J-JLEN-1+I )
756: 130 CONTINUE
757: ELSE
758: JLEN = MIN( KD, N-J )
759: DO 140 I = 1, JLEN
760: CSUMJ = CSUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
761: 140 CONTINUE
762: END IF
763: END IF
764: *
765: IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
766: *
767: * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
768: * was not used to scale the dotproduct.
769: *
770: X( J ) = X( J ) - CSUMJ
771: XJ = CABS1( X( J ) )
772: IF( NOUNIT ) THEN
773: *
774: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
775: *
776: TJJS = AB( MAIND, J )*TSCAL
777: ELSE
778: TJJS = TSCAL
779: IF( TSCAL.EQ.ONE )
780: $ GO TO 160
781: END IF
782: TJJ = CABS1( TJJS )
783: IF( TJJ.GT.SMLNUM ) THEN
784: *
785: * abs(A(j,j)) > SMLNUM:
786: *
787: IF( TJJ.LT.ONE ) THEN
788: IF( XJ.GT.TJJ*BIGNUM ) THEN
789: *
790: * Scale X by 1/abs(x(j)).
791: *
792: REC = ONE / XJ
793: CALL ZDSCAL( N, REC, X, 1 )
794: SCALE = SCALE*REC
795: XMAX = XMAX*REC
796: END IF
797: END IF
798: X( J ) = ZLADIV( X( J ), TJJS )
799: ELSE IF( TJJ.GT.ZERO ) THEN
800: *
801: * 0 < abs(A(j,j)) <= SMLNUM:
802: *
803: IF( XJ.GT.TJJ*BIGNUM ) THEN
804: *
805: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
806: *
807: REC = ( TJJ*BIGNUM ) / XJ
808: CALL ZDSCAL( N, REC, X, 1 )
809: SCALE = SCALE*REC
810: XMAX = XMAX*REC
811: END IF
812: X( J ) = ZLADIV( X( J ), TJJS )
813: ELSE
814: *
815: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
816: * scale = 0 and compute a solution to A**T *x = 0.
817: *
818: DO 150 I = 1, N
819: X( I ) = ZERO
820: 150 CONTINUE
821: X( J ) = ONE
822: SCALE = ZERO
823: XMAX = ZERO
824: END IF
825: 160 CONTINUE
826: ELSE
827: *
828: * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
829: * product has already been divided by 1/A(j,j).
830: *
831: X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
832: END IF
833: XMAX = MAX( XMAX, CABS1( X( J ) ) )
834: 170 CONTINUE
835: *
836: ELSE
837: *
838: * Solve A**H * x = b
839: *
840: DO 220 J = JFIRST, JLAST, JINC
841: *
842: * Compute x(j) = b(j) - sum A(k,j)*x(k).
843: * k<>j
844: *
845: XJ = CABS1( X( J ) )
846: USCAL = TSCAL
847: REC = ONE / MAX( XMAX, ONE )
848: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
849: *
850: * If x(j) could overflow, scale x by 1/(2*XMAX).
851: *
852: REC = REC*HALF
853: IF( NOUNIT ) THEN
854: TJJS = DCONJG( AB( MAIND, J ) )*TSCAL
855: ELSE
856: TJJS = TSCAL
857: END IF
858: TJJ = CABS1( TJJS )
859: IF( TJJ.GT.ONE ) THEN
860: *
861: * Divide by A(j,j) when scaling x if A(j,j) > 1.
862: *
863: REC = MIN( ONE, REC*TJJ )
864: USCAL = ZLADIV( USCAL, TJJS )
865: END IF
866: IF( REC.LT.ONE ) THEN
867: CALL ZDSCAL( N, REC, X, 1 )
868: SCALE = SCALE*REC
869: XMAX = XMAX*REC
870: END IF
871: END IF
872: *
873: CSUMJ = ZERO
874: IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
875: *
876: * If the scaling needed for A in the dot product is 1,
877: * call ZDOTC to perform the dot product.
878: *
879: IF( UPPER ) THEN
880: JLEN = MIN( KD, J-1 )
881: CSUMJ = ZDOTC( JLEN, AB( KD+1-JLEN, J ), 1,
882: $ X( J-JLEN ), 1 )
883: ELSE
884: JLEN = MIN( KD, N-J )
885: IF( JLEN.GT.1 )
886: $ CSUMJ = ZDOTC( JLEN, AB( 2, J ), 1, X( J+1 ),
887: $ 1 )
888: END IF
889: ELSE
890: *
891: * Otherwise, use in-line code for the dot product.
892: *
893: IF( UPPER ) THEN
894: JLEN = MIN( KD, J-1 )
895: DO 180 I = 1, JLEN
896: CSUMJ = CSUMJ + ( DCONJG( AB( KD+I-JLEN, J ) )*
897: $ USCAL )*X( J-JLEN-1+I )
898: 180 CONTINUE
899: ELSE
900: JLEN = MIN( KD, N-J )
901: DO 190 I = 1, JLEN
902: CSUMJ = CSUMJ + ( DCONJG( AB( I+1, J ) )*USCAL )
903: $ *X( J+I )
904: 190 CONTINUE
905: END IF
906: END IF
907: *
908: IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
909: *
910: * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
911: * was not used to scale the dotproduct.
912: *
913: X( J ) = X( J ) - CSUMJ
914: XJ = CABS1( X( J ) )
915: IF( NOUNIT ) THEN
916: *
917: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
918: *
919: TJJS = DCONJG( AB( MAIND, J ) )*TSCAL
920: ELSE
921: TJJS = TSCAL
922: IF( TSCAL.EQ.ONE )
923: $ GO TO 210
924: END IF
925: TJJ = CABS1( TJJS )
926: IF( TJJ.GT.SMLNUM ) THEN
927: *
928: * abs(A(j,j)) > SMLNUM:
929: *
930: IF( TJJ.LT.ONE ) THEN
931: IF( XJ.GT.TJJ*BIGNUM ) THEN
932: *
933: * Scale X by 1/abs(x(j)).
934: *
935: REC = ONE / XJ
936: CALL ZDSCAL( N, REC, X, 1 )
937: SCALE = SCALE*REC
938: XMAX = XMAX*REC
939: END IF
940: END IF
941: X( J ) = ZLADIV( X( J ), TJJS )
942: ELSE IF( TJJ.GT.ZERO ) THEN
943: *
944: * 0 < abs(A(j,j)) <= SMLNUM:
945: *
946: IF( XJ.GT.TJJ*BIGNUM ) THEN
947: *
948: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
949: *
950: REC = ( TJJ*BIGNUM ) / XJ
951: CALL ZDSCAL( N, REC, X, 1 )
952: SCALE = SCALE*REC
953: XMAX = XMAX*REC
954: END IF
955: X( J ) = ZLADIV( X( J ), TJJS )
956: ELSE
957: *
958: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
959: * scale = 0 and compute a solution to A**H *x = 0.
960: *
961: DO 200 I = 1, N
962: X( I ) = ZERO
963: 200 CONTINUE
964: X( J ) = ONE
965: SCALE = ZERO
966: XMAX = ZERO
967: END IF
968: 210 CONTINUE
969: ELSE
970: *
971: * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
972: * product has already been divided by 1/A(j,j).
973: *
974: X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
975: END IF
976: XMAX = MAX( XMAX, CABS1( X( J ) ) )
977: 220 CONTINUE
978: END IF
979: SCALE = SCALE / TSCAL
980: END IF
981: *
982: * Scale the column norms by 1/TSCAL for return.
983: *
984: IF( TSCAL.NE.ONE ) THEN
985: CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
986: END IF
987: *
988: RETURN
989: *
990: * End of ZLATBS
991: *
992: 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