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