001:       SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
002:      $                   CNORM, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          DIAG, NORMIN, TRANS, UPLO
011:       INTEGER            INFO, LDA, N
012:       DOUBLE PRECISION   SCALE
013: *     ..
014: *     .. Array Arguments ..
015:       DOUBLE PRECISION   A( LDA, * ), CNORM( * ), X( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DLATRS solves one of the triangular systems
022: *
023: *     A *x = s*b  or  A'*x = s*b
024: *
025: *  with scaling to prevent overflow.  Here A is an upper or lower
026: *  triangular matrix, A' denotes the transpose of A, x and b are
027: *  n-element vectors, and s is a scaling factor, usually less than
028: *  or equal to 1, chosen so that the components of x will be less than
029: *  the overflow threshold.  If the unscaled problem will not cause
030: *  overflow, the Level 2 BLAS routine DTRSV is called.  If the matrix A
031: *  is singular (A(j,j) = 0 for some j), then s is set to 0 and a
032: *  non-trivial solution to A*x = 0 is returned.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  UPLO    (input) CHARACTER*1
038: *          Specifies whether the matrix A is upper or lower triangular.
039: *          = 'U':  Upper triangular
040: *          = 'L':  Lower triangular
041: *
042: *  TRANS   (input) CHARACTER*1
043: *          Specifies the operation applied to A.
044: *          = 'N':  Solve A * x = s*b  (No transpose)
045: *          = 'T':  Solve A'* x = s*b  (Transpose)
046: *          = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)
047: *
048: *  DIAG    (input) CHARACTER*1
049: *          Specifies whether or not the matrix A is unit triangular.
050: *          = 'N':  Non-unit triangular
051: *          = 'U':  Unit triangular
052: *
053: *  NORMIN  (input) CHARACTER*1
054: *          Specifies whether CNORM has been set or not.
055: *          = 'Y':  CNORM contains the column norms on entry
056: *          = 'N':  CNORM is not set on entry.  On exit, the norms will
057: *                  be computed and stored in CNORM.
058: *
059: *  N       (input) INTEGER
060: *          The order of the matrix A.  N >= 0.
061: *
062: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
063: *          The triangular matrix A.  If UPLO = 'U', the leading n by n
064: *          upper triangular part of the array A contains the upper
065: *          triangular matrix, and the strictly lower triangular part of
066: *          A is not referenced.  If UPLO = 'L', the leading n by n lower
067: *          triangular part of the array A contains the lower triangular
068: *          matrix, and the strictly upper triangular part of A is not
069: *          referenced.  If DIAG = 'U', the diagonal elements of A are
070: *          also not referenced and are assumed to be 1.
071: *
072: *  LDA     (input) INTEGER
073: *          The leading dimension of the array A.  LDA >= max (1,N).
074: *
075: *  X       (input/output) DOUBLE PRECISION array, dimension (N)
076: *          On entry, the right hand side b of the triangular system.
077: *          On exit, X is overwritten by the solution vector x.
078: *
079: *  SCALE   (output) DOUBLE PRECISION
080: *          The scaling factor s for the triangular system
081: *             A * x = s*b  or  A'* x = s*b.
082: *          If SCALE = 0, the matrix A is singular or badly scaled, and
083: *          the vector x is an exact or approximate solution to A*x = 0.
084: *
085: *  CNORM   (input or output) DOUBLE PRECISION array, dimension (N)
086: *
087: *          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
088: *          contains the norm of the off-diagonal part of the j-th column
089: *          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
090: *          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
091: *          must be greater than or equal to the 1-norm.
092: *
093: *          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
094: *          returns the 1-norm of the offdiagonal part of the j-th column
095: *          of A.
096: *
097: *  INFO    (output) INTEGER
098: *          = 0:  successful exit
099: *          < 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
703: