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**T *x = s*b
24 *
25 * with scaling to prevent overflow. Here A is an upper or lower
26 * triangular matrix, A**T 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**T* x = s*b (Transpose)
46 * = 'C': Solve A**T* 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**T* 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**T*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]**T * 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]**T * 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**T * 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**T * 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**T*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
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**T *x = s*b
24 *
25 * with scaling to prevent overflow. Here A is an upper or lower
26 * triangular matrix, A**T 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**T* x = s*b (Transpose)
46 * = 'C': Solve A**T* 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**T* 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**T*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]**T * 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]**T * 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**T * 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**T * 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**T*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