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