Annotation of rpl/lapack/lapack/dlatrs.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLATRS
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLATRS + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlatrs.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlatrs.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlatrs.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLATRS( 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 A( LDA, * ), CNORM( * ), X( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DLATRS solves one of the triangular systems
! 40: *>
! 41: *> A *x = s*b or A**T *x = s*b
! 42: *>
! 43: *> with scaling to prevent overflow. Here A is an upper or lower
! 44: *> triangular matrix, A**T denotes the transpose of A, x and b are
! 45: *> n-element vectors, and s is a scaling factor, usually less than
! 46: *> or equal to 1, chosen so that the components of x will be less than
! 47: *> the overflow threshold. If the unscaled problem will not cause
! 48: *> overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A
! 49: *> is singular (A(j,j) = 0 for some j), then s is set to 0 and a
! 50: *> non-trivial solution to A*x = 0 is returned.
! 51: *> \endverbatim
! 52: *
! 53: * Arguments:
! 54: * ==========
! 55: *
! 56: *> \param[in] UPLO
! 57: *> \verbatim
! 58: *> UPLO is CHARACTER*1
! 59: *> Specifies whether the matrix A is upper or lower triangular.
! 60: *> = 'U': Upper triangular
! 61: *> = 'L': Lower triangular
! 62: *> \endverbatim
! 63: *>
! 64: *> \param[in] TRANS
! 65: *> \verbatim
! 66: *> TRANS is CHARACTER*1
! 67: *> Specifies the operation applied to A.
! 68: *> = 'N': Solve A * x = s*b (No transpose)
! 69: *> = 'T': Solve A**T* x = s*b (Transpose)
! 70: *> = 'C': Solve A**T* x = s*b (Conjugate transpose = Transpose)
! 71: *> \endverbatim
! 72: *>
! 73: *> \param[in] DIAG
! 74: *> \verbatim
! 75: *> DIAG is CHARACTER*1
! 76: *> Specifies whether or not the matrix A is unit triangular.
! 77: *> = 'N': Non-unit triangular
! 78: *> = 'U': Unit triangular
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] NORMIN
! 82: *> \verbatim
! 83: *> NORMIN is CHARACTER*1
! 84: *> Specifies whether CNORM has been set or not.
! 85: *> = 'Y': CNORM contains the column norms on entry
! 86: *> = 'N': CNORM is not set on entry. On exit, the norms will
! 87: *> be computed and stored in CNORM.
! 88: *> \endverbatim
! 89: *>
! 90: *> \param[in] N
! 91: *> \verbatim
! 92: *> N is INTEGER
! 93: *> The order of the matrix A. N >= 0.
! 94: *> \endverbatim
! 95: *>
! 96: *> \param[in] A
! 97: *> \verbatim
! 98: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 99: *> The triangular matrix A. If UPLO = 'U', the leading n by n
! 100: *> upper triangular part of the array A contains the upper
! 101: *> triangular matrix, and the strictly lower triangular part of
! 102: *> A is not referenced. If UPLO = 'L', the leading n by n lower
! 103: *> triangular part of the array A contains the lower triangular
! 104: *> matrix, and the strictly upper triangular part of A is not
! 105: *> referenced. If DIAG = 'U', the diagonal elements of A are
! 106: *> also not referenced and are assumed to be 1.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[in] LDA
! 110: *> \verbatim
! 111: *> LDA is INTEGER
! 112: *> The leading dimension of the array A. LDA >= max (1,N).
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in,out] X
! 116: *> \verbatim
! 117: *> X is DOUBLE PRECISION array, dimension (N)
! 118: *> On entry, the right hand side b of the triangular system.
! 119: *> On exit, X is overwritten by the solution vector x.
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[out] SCALE
! 123: *> \verbatim
! 124: *> SCALE is DOUBLE PRECISION
! 125: *> The scaling factor s for the triangular system
! 126: *> A * x = s*b or A**T* x = s*b.
! 127: *> If SCALE = 0, the matrix A is singular or badly scaled, and
! 128: *> the vector x is an exact or approximate solution to A*x = 0.
! 129: *> \endverbatim
! 130: *>
! 131: *> \param[in,out] CNORM
! 132: *> \verbatim
! 133: *> CNORM is or output) DOUBLE PRECISION array, dimension (N)
! 134: *>
! 135: *> If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
! 136: *> contains the norm of the off-diagonal part of the j-th column
! 137: *> of A. If TRANS = 'N', CNORM(j) must be greater than or equal
! 138: *> to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
! 139: *> must be greater than or equal to the 1-norm.
! 140: *>
! 141: *> If NORMIN = 'N', CNORM is an output argument and CNORM(j)
! 142: *> returns the 1-norm of the offdiagonal part of the j-th column
! 143: *> of A.
! 144: *> \endverbatim
! 145: *>
! 146: *> \param[out] INFO
! 147: *> \verbatim
! 148: *> INFO is INTEGER
! 149: *> = 0: successful exit
! 150: *> < 0: if INFO = -k, the k-th argument had an illegal value
! 151: *> \endverbatim
! 152: *
! 153: * Authors:
! 154: * ========
! 155: *
! 156: *> \author Univ. of Tennessee
! 157: *> \author Univ. of California Berkeley
! 158: *> \author Univ. of Colorado Denver
! 159: *> \author NAG Ltd.
! 160: *
! 161: *> \date November 2011
! 162: *
! 163: *> \ingroup doubleOTHERauxiliary
! 164: *
! 165: *> \par Further Details:
! 166: * =====================
! 167: *>
! 168: *> \verbatim
! 169: *>
! 170: *> A rough bound on x is computed; if that is less than overflow, DTRSV
! 171: *> is called, otherwise, specific code is used which checks for possible
! 172: *> overflow or divide-by-zero at every operation.
! 173: *>
! 174: *> A columnwise scheme is used for solving A*x = b. The basic algorithm
! 175: *> if A is lower triangular is
! 176: *>
! 177: *> x[1:n] := b[1:n]
! 178: *> for j = 1, ..., n
! 179: *> x(j) := x(j) / A(j,j)
! 180: *> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
! 181: *> end
! 182: *>
! 183: *> Define bounds on the components of x after j iterations of the loop:
! 184: *> M(j) = bound on x[1:j]
! 185: *> G(j) = bound on x[j+1:n]
! 186: *> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
! 187: *>
! 188: *> Then for iteration j+1 we have
! 189: *> M(j+1) <= G(j) / | A(j+1,j+1) |
! 190: *> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
! 191: *> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
! 192: *>
! 193: *> where CNORM(j+1) is greater than or equal to the infinity-norm of
! 194: *> column j+1 of A, not counting the diagonal. Hence
! 195: *>
! 196: *> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
! 197: *> 1<=i<=j
! 198: *> and
! 199: *>
! 200: *> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
! 201: *> 1<=i< j
! 202: *>
! 203: *> Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the
! 204: *> reciprocal of the largest M(j), j=1,..,n, is larger than
! 205: *> max(underflow, 1/overflow).
! 206: *>
! 207: *> The bound on x(j) is also used to determine when a step in the
! 208: *> columnwise method can be performed without fear of overflow. If
! 209: *> the computed bound is greater than a large constant, x is scaled to
! 210: *> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
! 211: *> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
! 212: *>
! 213: *> Similarly, a row-wise scheme is used to solve A**T*x = b. The basic
! 214: *> algorithm for A upper triangular is
! 215: *>
! 216: *> for j = 1, ..., n
! 217: *> x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
! 218: *> end
! 219: *>
! 220: *> We simultaneously compute two bounds
! 221: *> G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
! 222: *> M(j) = bound on x(i), 1<=i<=j
! 223: *>
! 224: *> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
! 225: *> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
! 226: *> Then the bound on x(j) is
! 227: *>
! 228: *> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
! 229: *>
! 230: *> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
! 231: *> 1<=i<=j
! 232: *>
! 233: *> and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater
! 234: *> than max(underflow, 1/overflow).
! 235: *> \endverbatim
! 236: *>
! 237: * =====================================================================
1.1 bertrand 238: SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
239: $ CNORM, INFO )
240: *
1.9 ! bertrand 241: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 242: * -- LAPACK is a software package provided by Univ. of Tennessee, --
243: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 244: * November 2011
1.1 bertrand 245: *
246: * .. Scalar Arguments ..
247: CHARACTER DIAG, NORMIN, TRANS, UPLO
248: INTEGER INFO, LDA, N
249: DOUBLE PRECISION SCALE
250: * ..
251: * .. Array Arguments ..
252: DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * )
253: * ..
254: *
255: * =====================================================================
256: *
257: * .. Parameters ..
258: DOUBLE PRECISION ZERO, HALF, ONE
259: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.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, SUMJ, TJJ, TJJS,
265: $ TMAX, TSCAL, USCAL, XBND, XJ, XMAX
266: * ..
267: * .. External Functions ..
268: LOGICAL LSAME
269: INTEGER IDAMAX
270: DOUBLE PRECISION DASUM, DDOT, DLAMCH
271: EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
272: * ..
273: * .. External Subroutines ..
274: EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA
275: * ..
276: * .. Intrinsic Functions ..
277: INTRINSIC ABS, MAX, MIN
278: * ..
279: * .. Executable Statements ..
280: *
281: INFO = 0
282: UPPER = LSAME( UPLO, 'U' )
283: NOTRAN = LSAME( TRANS, 'N' )
284: NOUNIT = LSAME( DIAG, 'N' )
285: *
286: * Test the input parameters.
287: *
288: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
289: INFO = -1
290: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
291: $ LSAME( TRANS, 'C' ) ) THEN
292: INFO = -2
293: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
294: INFO = -3
295: ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
296: $ LSAME( NORMIN, 'N' ) ) THEN
297: INFO = -4
298: ELSE IF( N.LT.0 ) THEN
299: INFO = -5
300: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
301: INFO = -7
302: END IF
303: IF( INFO.NE.0 ) THEN
304: CALL XERBLA( 'DLATRS', -INFO )
305: RETURN
306: END IF
307: *
308: * Quick return if possible
309: *
310: IF( N.EQ.0 )
311: $ RETURN
312: *
313: * Determine machine dependent parameters to control overflow.
314: *
315: SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
316: BIGNUM = ONE / SMLNUM
317: SCALE = ONE
318: *
319: IF( LSAME( NORMIN, 'N' ) ) THEN
320: *
321: * Compute the 1-norm of each column, not including the diagonal.
322: *
323: IF( UPPER ) THEN
324: *
325: * A is upper triangular.
326: *
327: DO 10 J = 1, N
328: CNORM( J ) = DASUM( J-1, A( 1, J ), 1 )
329: 10 CONTINUE
330: ELSE
331: *
332: * A is lower triangular.
333: *
334: DO 20 J = 1, N - 1
335: CNORM( J ) = DASUM( N-J, A( J+1, J ), 1 )
336: 20 CONTINUE
337: CNORM( N ) = ZERO
338: END IF
339: END IF
340: *
341: * Scale the column norms by TSCAL if the maximum element in CNORM is
342: * greater than BIGNUM.
343: *
344: IMAX = IDAMAX( N, CNORM, 1 )
345: TMAX = CNORM( IMAX )
346: IF( TMAX.LE.BIGNUM ) THEN
347: TSCAL = ONE
348: ELSE
349: TSCAL = ONE / ( SMLNUM*TMAX )
350: CALL DSCAL( N, TSCAL, CNORM, 1 )
351: END IF
352: *
353: * Compute a bound on the computed solution vector to see if the
354: * Level 2 BLAS routine DTRSV can be used.
355: *
356: J = IDAMAX( N, X, 1 )
357: XMAX = ABS( X( J ) )
358: XBND = XMAX
359: IF( NOTRAN ) THEN
360: *
361: * Compute the growth in A * x = b.
362: *
363: IF( UPPER ) THEN
364: JFIRST = N
365: JLAST = 1
366: JINC = -1
367: ELSE
368: JFIRST = 1
369: JLAST = N
370: JINC = 1
371: END IF
372: *
373: IF( TSCAL.NE.ONE ) THEN
374: GROW = ZERO
375: GO TO 50
376: END IF
377: *
378: IF( NOUNIT ) THEN
379: *
380: * A is non-unit triangular.
381: *
382: * Compute GROW = 1/G(j) and XBND = 1/M(j).
383: * Initially, G(0) = max{x(i), i=1,...,n}.
384: *
385: GROW = ONE / MAX( XBND, SMLNUM )
386: XBND = GROW
387: DO 30 J = JFIRST, JLAST, JINC
388: *
389: * Exit the loop if the growth factor is too small.
390: *
391: IF( GROW.LE.SMLNUM )
392: $ GO TO 50
393: *
394: * M(j) = G(j-1) / abs(A(j,j))
395: *
396: TJJ = ABS( A( J, J ) )
397: XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
398: IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
399: *
400: * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
401: *
402: GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
403: ELSE
404: *
405: * G(j) could overflow, set GROW to 0.
406: *
407: GROW = ZERO
408: END IF
409: 30 CONTINUE
410: GROW = XBND
411: ELSE
412: *
413: * A is unit triangular.
414: *
415: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
416: *
417: GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
418: DO 40 J = JFIRST, JLAST, JINC
419: *
420: * Exit the loop if the growth factor is too small.
421: *
422: IF( GROW.LE.SMLNUM )
423: $ GO TO 50
424: *
425: * G(j) = G(j-1)*( 1 + CNORM(j) )
426: *
427: GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
428: 40 CONTINUE
429: END IF
430: 50 CONTINUE
431: *
432: ELSE
433: *
1.8 bertrand 434: * Compute the growth in A**T * x = b.
1.1 bertrand 435: *
436: IF( UPPER ) THEN
437: JFIRST = 1
438: JLAST = N
439: JINC = 1
440: ELSE
441: JFIRST = N
442: JLAST = 1
443: JINC = -1
444: END IF
445: *
446: IF( TSCAL.NE.ONE ) THEN
447: GROW = ZERO
448: GO TO 80
449: END IF
450: *
451: IF( NOUNIT ) THEN
452: *
453: * A is non-unit triangular.
454: *
455: * Compute GROW = 1/G(j) and XBND = 1/M(j).
456: * Initially, M(0) = max{x(i), i=1,...,n}.
457: *
458: GROW = ONE / MAX( XBND, SMLNUM )
459: XBND = GROW
460: DO 60 J = JFIRST, JLAST, JINC
461: *
462: * Exit the loop if the growth factor is too small.
463: *
464: IF( GROW.LE.SMLNUM )
465: $ GO TO 80
466: *
467: * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
468: *
469: XJ = ONE + CNORM( J )
470: GROW = MIN( GROW, XBND / XJ )
471: *
472: * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
473: *
474: TJJ = ABS( A( J, J ) )
475: IF( XJ.GT.TJJ )
476: $ XBND = XBND*( TJJ / XJ )
477: 60 CONTINUE
478: GROW = MIN( GROW, XBND )
479: ELSE
480: *
481: * A is unit triangular.
482: *
483: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
484: *
485: GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
486: DO 70 J = JFIRST, JLAST, JINC
487: *
488: * Exit the loop if the growth factor is too small.
489: *
490: IF( GROW.LE.SMLNUM )
491: $ GO TO 80
492: *
493: * G(j) = ( 1 + CNORM(j) )*G(j-1)
494: *
495: XJ = ONE + CNORM( J )
496: GROW = GROW / XJ
497: 70 CONTINUE
498: END IF
499: 80 CONTINUE
500: END IF
501: *
502: IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
503: *
504: * Use the Level 2 BLAS solve if the reciprocal of the bound on
505: * elements of X is not too small.
506: *
507: CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
508: ELSE
509: *
510: * Use a Level 1 BLAS solve, scaling intermediate results.
511: *
512: IF( XMAX.GT.BIGNUM ) THEN
513: *
514: * Scale X so that its components are less than or equal to
515: * BIGNUM in absolute value.
516: *
517: SCALE = BIGNUM / XMAX
518: CALL DSCAL( N, SCALE, X, 1 )
519: XMAX = BIGNUM
520: END IF
521: *
522: IF( NOTRAN ) THEN
523: *
524: * Solve A * x = b
525: *
526: DO 110 J = JFIRST, JLAST, JINC
527: *
528: * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
529: *
530: XJ = ABS( X( J ) )
531: IF( NOUNIT ) THEN
532: TJJS = A( J, J )*TSCAL
533: ELSE
534: TJJS = TSCAL
535: IF( TSCAL.EQ.ONE )
536: $ GO TO 100
537: END IF
538: TJJ = ABS( TJJS )
539: IF( TJJ.GT.SMLNUM ) THEN
540: *
541: * abs(A(j,j)) > SMLNUM:
542: *
543: IF( TJJ.LT.ONE ) THEN
544: IF( XJ.GT.TJJ*BIGNUM ) THEN
545: *
546: * Scale x by 1/b(j).
547: *
548: REC = ONE / XJ
549: CALL DSCAL( N, REC, X, 1 )
550: SCALE = SCALE*REC
551: XMAX = XMAX*REC
552: END IF
553: END IF
554: X( J ) = X( J ) / TJJS
555: XJ = ABS( X( J ) )
556: ELSE IF( TJJ.GT.ZERO ) THEN
557: *
558: * 0 < abs(A(j,j)) <= SMLNUM:
559: *
560: IF( XJ.GT.TJJ*BIGNUM ) THEN
561: *
562: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
563: * to avoid overflow when dividing by A(j,j).
564: *
565: REC = ( TJJ*BIGNUM ) / XJ
566: IF( CNORM( J ).GT.ONE ) THEN
567: *
568: * Scale by 1/CNORM(j) to avoid overflow when
569: * multiplying x(j) times column j.
570: *
571: REC = REC / CNORM( J )
572: END IF
573: CALL DSCAL( N, REC, X, 1 )
574: SCALE = SCALE*REC
575: XMAX = XMAX*REC
576: END IF
577: X( J ) = X( J ) / TJJS
578: XJ = ABS( X( J ) )
579: ELSE
580: *
581: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
582: * scale = 0, and compute a solution to A*x = 0.
583: *
584: DO 90 I = 1, N
585: X( I ) = ZERO
586: 90 CONTINUE
587: X( J ) = ONE
588: XJ = ONE
589: SCALE = ZERO
590: XMAX = ZERO
591: END IF
592: 100 CONTINUE
593: *
594: * Scale x if necessary to avoid overflow when adding a
595: * multiple of column j of A.
596: *
597: IF( XJ.GT.ONE ) THEN
598: REC = ONE / XJ
599: IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
600: *
601: * Scale x by 1/(2*abs(x(j))).
602: *
603: REC = REC*HALF
604: CALL DSCAL( N, REC, X, 1 )
605: SCALE = SCALE*REC
606: END IF
607: ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
608: *
609: * Scale x by 1/2.
610: *
611: CALL DSCAL( N, HALF, X, 1 )
612: SCALE = SCALE*HALF
613: END IF
614: *
615: IF( UPPER ) THEN
616: IF( J.GT.1 ) THEN
617: *
618: * Compute the update
619: * x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
620: *
621: CALL DAXPY( J-1, -X( J )*TSCAL, A( 1, J ), 1, X,
622: $ 1 )
623: I = IDAMAX( J-1, X, 1 )
624: XMAX = ABS( X( I ) )
625: END IF
626: ELSE
627: IF( J.LT.N ) THEN
628: *
629: * Compute the update
630: * x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
631: *
632: CALL DAXPY( N-J, -X( J )*TSCAL, A( J+1, J ), 1,
633: $ X( J+1 ), 1 )
634: I = J + IDAMAX( N-J, X( J+1 ), 1 )
635: XMAX = ABS( X( I ) )
636: END IF
637: END IF
638: 110 CONTINUE
639: *
640: ELSE
641: *
1.8 bertrand 642: * Solve A**T * x = b
1.1 bertrand 643: *
644: DO 160 J = JFIRST, JLAST, JINC
645: *
646: * Compute x(j) = b(j) - sum A(k,j)*x(k).
647: * k<>j
648: *
649: XJ = ABS( X( J ) )
650: USCAL = TSCAL
651: REC = ONE / MAX( XMAX, ONE )
652: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
653: *
654: * If x(j) could overflow, scale x by 1/(2*XMAX).
655: *
656: REC = REC*HALF
657: IF( NOUNIT ) THEN
658: TJJS = A( J, J )*TSCAL
659: ELSE
660: TJJS = TSCAL
661: END IF
662: TJJ = ABS( TJJS )
663: IF( TJJ.GT.ONE ) THEN
664: *
665: * Divide by A(j,j) when scaling x if A(j,j) > 1.
666: *
667: REC = MIN( ONE, REC*TJJ )
668: USCAL = USCAL / TJJS
669: END IF
670: IF( REC.LT.ONE ) THEN
671: CALL DSCAL( N, REC, X, 1 )
672: SCALE = SCALE*REC
673: XMAX = XMAX*REC
674: END IF
675: END IF
676: *
677: SUMJ = ZERO
678: IF( USCAL.EQ.ONE ) THEN
679: *
680: * If the scaling needed for A in the dot product is 1,
681: * call DDOT to perform the dot product.
682: *
683: IF( UPPER ) THEN
684: SUMJ = DDOT( J-1, A( 1, J ), 1, X, 1 )
685: ELSE IF( J.LT.N ) THEN
686: SUMJ = DDOT( N-J, A( J+1, J ), 1, X( J+1 ), 1 )
687: END IF
688: ELSE
689: *
690: * Otherwise, use in-line code for the dot product.
691: *
692: IF( UPPER ) THEN
693: DO 120 I = 1, J - 1
694: SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I )
695: 120 CONTINUE
696: ELSE IF( J.LT.N ) THEN
697: DO 130 I = J + 1, N
698: SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I )
699: 130 CONTINUE
700: END IF
701: END IF
702: *
703: IF( USCAL.EQ.TSCAL ) THEN
704: *
705: * Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
706: * was not used to scale the dotproduct.
707: *
708: X( J ) = X( J ) - SUMJ
709: XJ = ABS( X( J ) )
710: IF( NOUNIT ) THEN
711: TJJS = A( J, J )*TSCAL
712: ELSE
713: TJJS = TSCAL
714: IF( TSCAL.EQ.ONE )
715: $ GO TO 150
716: END IF
717: *
718: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
719: *
720: TJJ = ABS( TJJS )
721: IF( TJJ.GT.SMLNUM ) THEN
722: *
723: * abs(A(j,j)) > SMLNUM:
724: *
725: IF( TJJ.LT.ONE ) THEN
726: IF( XJ.GT.TJJ*BIGNUM ) THEN
727: *
728: * Scale X by 1/abs(x(j)).
729: *
730: REC = ONE / XJ
731: CALL DSCAL( N, REC, X, 1 )
732: SCALE = SCALE*REC
733: XMAX = XMAX*REC
734: END IF
735: END IF
736: X( J ) = X( J ) / TJJS
737: ELSE IF( TJJ.GT.ZERO ) THEN
738: *
739: * 0 < abs(A(j,j)) <= SMLNUM:
740: *
741: IF( XJ.GT.TJJ*BIGNUM ) THEN
742: *
743: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
744: *
745: REC = ( TJJ*BIGNUM ) / XJ
746: CALL DSCAL( N, REC, X, 1 )
747: SCALE = SCALE*REC
748: XMAX = XMAX*REC
749: END IF
750: X( J ) = X( J ) / TJJS
751: ELSE
752: *
753: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
1.8 bertrand 754: * scale = 0, and compute a solution to A**T*x = 0.
1.1 bertrand 755: *
756: DO 140 I = 1, N
757: X( I ) = ZERO
758: 140 CONTINUE
759: X( J ) = ONE
760: SCALE = ZERO
761: XMAX = ZERO
762: END IF
763: 150 CONTINUE
764: ELSE
765: *
766: * Compute x(j) := x(j) / A(j,j) - sumj if the dot
767: * product has already been divided by 1/A(j,j).
768: *
769: X( J ) = X( J ) / TJJS - SUMJ
770: END IF
771: XMAX = MAX( XMAX, ABS( X( J ) ) )
772: 160 CONTINUE
773: END IF
774: SCALE = SCALE / TSCAL
775: END IF
776: *
777: * Scale the column norms by 1/TSCAL for return.
778: *
779: IF( TSCAL.NE.ONE ) THEN
780: CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
781: END IF
782: *
783: RETURN
784: *
785: * End of DLATRS
786: *
787: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlatrs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack