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