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