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