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