```001:       SUBROUTINE CLATRS( 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:       REAL               SCALE
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               CNORM( * )
015:       COMPLEX            A( LDA, * ), X( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CLATRS solves one of the triangular systems
022: *
023: *     A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,
024: *
025: *  with scaling to prevent overflow.  Here A is an upper or lower
026: *  triangular matrix, A**T denotes the transpose of A, A**H denotes the
027: *  conjugate transpose of A, x and b are n-element vectors, and s is a
028: *  scaling factor, usually less than or equal to 1, chosen so that the
029: *  components of x will be less than the overflow threshold.  If the
030: *  unscaled problem will not cause overflow, the Level 2 BLAS routine
031: *  CTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
032: *  then s is set to 0 and a 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**T * x = s*b  (Transpose)
046: *          = 'C':  Solve A**H * x = s*b  (Conjugate 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) COMPLEX 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) COMPLEX 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) REAL
080: *          The scaling factor s for the triangular system
081: *             A * x = s*b,  A**T * x = s*b,  or  A**H * 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) REAL 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, CTRSV
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 CTRSV 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  or
148: *  A**H *x = b.  The basic 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 CTRSV if 1/M(n) and 1/G(n) are both greater
168: *  than max(underflow, 1/overflow).
169: *
170: *  =====================================================================
171: *
172: *     .. Parameters ..
173:       REAL               ZERO, HALF, ONE, TWO
174:       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0,
175:      \$                   TWO = 2.0E+0 )
176: *     ..
177: *     .. Local Scalars ..
178:       LOGICAL            NOTRAN, NOUNIT, UPPER
179:       INTEGER            I, IMAX, J, JFIRST, JINC, JLAST
180:       REAL               BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
181:      \$                   XBND, XJ, XMAX
182:       COMPLEX            CSUMJ, TJJS, USCAL, ZDUM
183: *     ..
184: *     .. External Functions ..
185:       LOGICAL            LSAME
186:       INTEGER            ICAMAX, ISAMAX
187:       REAL               SCASUM, SLAMCH
189:       EXTERNAL           LSAME, ICAMAX, ISAMAX, SCASUM, SLAMCH, CDOTC,
191: *     ..
192: *     .. External Subroutines ..
193:       EXTERNAL           CAXPY, CSSCAL, CTRSV, SLABAD, SSCAL, XERBLA
194: *     ..
195: *     .. Intrinsic Functions ..
196:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
197: *     ..
198: *     .. Statement Functions ..
199:       REAL               CABS1, CABS2
200: *     ..
201: *     .. Statement Function definitions ..
202:       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
203:       CABS2( ZDUM ) = ABS( REAL( ZDUM ) / 2. ) +
204:      \$                ABS( AIMAG( ZDUM ) / 2. )
205: *     ..
206: *     .. Executable Statements ..
207: *
208:       INFO = 0
209:       UPPER = LSAME( UPLO, 'U' )
210:       NOTRAN = LSAME( TRANS, 'N' )
211:       NOUNIT = LSAME( DIAG, 'N' )
212: *
213: *     Test the input parameters.
214: *
215:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
216:          INFO = -1
217:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
218:      \$         LSAME( TRANS, 'C' ) ) THEN
219:          INFO = -2
220:       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
221:          INFO = -3
222:       ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
223:      \$         LSAME( NORMIN, 'N' ) ) THEN
224:          INFO = -4
225:       ELSE IF( N.LT.0 ) THEN
226:          INFO = -5
227:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
228:          INFO = -7
229:       END IF
230:       IF( INFO.NE.0 ) THEN
231:          CALL XERBLA( 'CLATRS', -INFO )
232:          RETURN
233:       END IF
234: *
235: *     Quick return if possible
236: *
237:       IF( N.EQ.0 )
238:      \$   RETURN
239: *
240: *     Determine machine dependent parameters to control overflow.
241: *
242:       SMLNUM = SLAMCH( 'Safe minimum' )
243:       BIGNUM = ONE / SMLNUM
244:       CALL SLABAD( SMLNUM, BIGNUM )
245:       SMLNUM = SMLNUM / SLAMCH( 'Precision' )
246:       BIGNUM = ONE / SMLNUM
247:       SCALE = ONE
248: *
249:       IF( LSAME( NORMIN, 'N' ) ) THEN
250: *
251: *        Compute the 1-norm of each column, not including the diagonal.
252: *
253:          IF( UPPER ) THEN
254: *
255: *           A is upper triangular.
256: *
257:             DO 10 J = 1, N
258:                CNORM( J ) = SCASUM( J-1, A( 1, J ), 1 )
259:    10       CONTINUE
260:          ELSE
261: *
262: *           A is lower triangular.
263: *
264:             DO 20 J = 1, N - 1
265:                CNORM( J ) = SCASUM( N-J, A( J+1, J ), 1 )
266:    20       CONTINUE
267:             CNORM( N ) = ZERO
268:          END IF
269:       END IF
270: *
271: *     Scale the column norms by TSCAL if the maximum element in CNORM is
272: *     greater than BIGNUM/2.
273: *
274:       IMAX = ISAMAX( N, CNORM, 1 )
275:       TMAX = CNORM( IMAX )
276:       IF( TMAX.LE.BIGNUM*HALF ) THEN
277:          TSCAL = ONE
278:       ELSE
279:          TSCAL = HALF / ( SMLNUM*TMAX )
280:          CALL SSCAL( N, TSCAL, CNORM, 1 )
281:       END IF
282: *
283: *     Compute a bound on the computed solution vector to see if the
284: *     Level 2 BLAS routine CTRSV can be used.
285: *
286:       XMAX = ZERO
287:       DO 30 J = 1, N
288:          XMAX = MAX( XMAX, CABS2( X( J ) ) )
289:    30 CONTINUE
290:       XBND = XMAX
291: *
292:       IF( NOTRAN ) THEN
293: *
294: *        Compute the growth in A * x = b.
295: *
296:          IF( UPPER ) THEN
297:             JFIRST = N
298:             JLAST = 1
299:             JINC = -1
300:          ELSE
301:             JFIRST = 1
302:             JLAST = N
303:             JINC = 1
304:          END IF
305: *
306:          IF( TSCAL.NE.ONE ) THEN
307:             GROW = ZERO
308:             GO TO 60
309:          END IF
310: *
311:          IF( NOUNIT ) THEN
312: *
313: *           A is non-unit triangular.
314: *
315: *           Compute GROW = 1/G(j) and XBND = 1/M(j).
316: *           Initially, G(0) = max{x(i), i=1,...,n}.
317: *
318:             GROW = HALF / MAX( XBND, SMLNUM )
319:             XBND = GROW
320:             DO 40 J = JFIRST, JLAST, JINC
321: *
322: *              Exit the loop if the growth factor is too small.
323: *
324:                IF( GROW.LE.SMLNUM )
325:      \$            GO TO 60
326: *
327:                TJJS = A( J, J )
328:                TJJ = CABS1( TJJS )
329: *
330:                IF( TJJ.GE.SMLNUM ) THEN
331: *
332: *                 M(j) = G(j-1) / abs(A(j,j))
333: *
334:                   XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
335:                ELSE
336: *
337: *                 M(j) could overflow, set XBND to 0.
338: *
339:                   XBND = ZERO
340:                END IF
341: *
342:                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
343: *
344: *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
345: *
346:                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
347:                ELSE
348: *
349: *                 G(j) could overflow, set GROW to 0.
350: *
351:                   GROW = ZERO
352:                END IF
353:    40       CONTINUE
354:             GROW = XBND
355:          ELSE
356: *
357: *           A is unit triangular.
358: *
359: *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
360: *
361:             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
362:             DO 50 J = JFIRST, JLAST, JINC
363: *
364: *              Exit the loop if the growth factor is too small.
365: *
366:                IF( GROW.LE.SMLNUM )
367:      \$            GO TO 60
368: *
369: *              G(j) = G(j-1)*( 1 + CNORM(j) )
370: *
371:                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
372:    50       CONTINUE
373:          END IF
374:    60    CONTINUE
375: *
376:       ELSE
377: *
378: *        Compute the growth in A**T * x = b  or  A**H * x = b.
379: *
380:          IF( UPPER ) THEN
381:             JFIRST = 1
382:             JLAST = N
383:             JINC = 1
384:          ELSE
385:             JFIRST = N
386:             JLAST = 1
387:             JINC = -1
388:          END IF
389: *
390:          IF( TSCAL.NE.ONE ) THEN
391:             GROW = ZERO
392:             GO TO 90
393:          END IF
394: *
395:          IF( NOUNIT ) THEN
396: *
397: *           A is non-unit triangular.
398: *
399: *           Compute GROW = 1/G(j) and XBND = 1/M(j).
400: *           Initially, M(0) = max{x(i), i=1,...,n}.
401: *
402:             GROW = HALF / MAX( XBND, SMLNUM )
403:             XBND = GROW
404:             DO 70 J = JFIRST, JLAST, JINC
405: *
406: *              Exit the loop if the growth factor is too small.
407: *
408:                IF( GROW.LE.SMLNUM )
409:      \$            GO TO 90
410: *
411: *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
412: *
413:                XJ = ONE + CNORM( J )
414:                GROW = MIN( GROW, XBND / XJ )
415: *
416:                TJJS = A( J, J )
417:                TJJ = CABS1( TJJS )
418: *
419:                IF( TJJ.GE.SMLNUM ) THEN
420: *
421: *                 M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
422: *
423:                   IF( XJ.GT.TJJ )
424:      \$               XBND = XBND*( TJJ / XJ )
425:                ELSE
426: *
427: *                 M(j) could overflow, set XBND to 0.
428: *
429:                   XBND = ZERO
430:                END IF
431:    70       CONTINUE
432:             GROW = MIN( GROW, XBND )
433:          ELSE
434: *
435: *           A is unit triangular.
436: *
437: *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
438: *
439:             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
440:             DO 80 J = JFIRST, JLAST, JINC
441: *
442: *              Exit the loop if the growth factor is too small.
443: *
444:                IF( GROW.LE.SMLNUM )
445:      \$            GO TO 90
446: *
447: *              G(j) = ( 1 + CNORM(j) )*G(j-1)
448: *
449:                XJ = ONE + CNORM( J )
450:                GROW = GROW / XJ
451:    80       CONTINUE
452:          END IF
453:    90    CONTINUE
454:       END IF
455: *
456:       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
457: *
458: *        Use the Level 2 BLAS solve if the reciprocal of the bound on
459: *        elements of X is not too small.
460: *
461:          CALL CTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
462:       ELSE
463: *
464: *        Use a Level 1 BLAS solve, scaling intermediate results.
465: *
466:          IF( XMAX.GT.BIGNUM*HALF ) THEN
467: *
468: *           Scale X so that its components are less than or equal to
469: *           BIGNUM in absolute value.
470: *
471:             SCALE = ( BIGNUM*HALF ) / XMAX
472:             CALL CSSCAL( N, SCALE, X, 1 )
473:             XMAX = BIGNUM
474:          ELSE
475:             XMAX = XMAX*TWO
476:          END IF
477: *
478:          IF( NOTRAN ) THEN
479: *
480: *           Solve A * x = b
481: *
482:             DO 110 J = JFIRST, JLAST, JINC
483: *
484: *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
485: *
486:                XJ = CABS1( X( J ) )
487:                IF( NOUNIT ) THEN
488:                   TJJS = A( J, J )*TSCAL
489:                ELSE
490:                   TJJS = TSCAL
491:                   IF( TSCAL.EQ.ONE )
492:      \$               GO TO 105
493:                END IF
494:                   TJJ = CABS1( TJJS )
495:                   IF( TJJ.GT.SMLNUM ) THEN
496: *
497: *                    abs(A(j,j)) > SMLNUM:
498: *
499:                      IF( TJJ.LT.ONE ) THEN
500:                         IF( XJ.GT.TJJ*BIGNUM ) THEN
501: *
502: *                          Scale x by 1/b(j).
503: *
504:                            REC = ONE / XJ
505:                            CALL CSSCAL( N, REC, X, 1 )
506:                            SCALE = SCALE*REC
507:                            XMAX = XMAX*REC
508:                         END IF
509:                      END IF
510:                      X( J ) = CLADIV( X( J ), TJJS )
511:                      XJ = CABS1( X( J ) )
512:                   ELSE IF( TJJ.GT.ZERO ) THEN
513: *
514: *                    0 < abs(A(j,j)) <= SMLNUM:
515: *
516:                      IF( XJ.GT.TJJ*BIGNUM ) THEN
517: *
518: *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
519: *                       to avoid overflow when dividing by A(j,j).
520: *
521:                         REC = ( TJJ*BIGNUM ) / XJ
522:                         IF( CNORM( J ).GT.ONE ) THEN
523: *
524: *                          Scale by 1/CNORM(j) to avoid overflow when
525: *                          multiplying x(j) times column j.
526: *
527:                            REC = REC / CNORM( J )
528:                         END IF
529:                         CALL CSSCAL( N, REC, X, 1 )
530:                         SCALE = SCALE*REC
531:                         XMAX = XMAX*REC
532:                      END IF
533:                      X( J ) = CLADIV( X( J ), TJJS )
534:                      XJ = CABS1( X( J ) )
535:                   ELSE
536: *
537: *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
538: *                    scale = 0, and compute a solution to A*x = 0.
539: *
540:                      DO 100 I = 1, N
541:                         X( I ) = ZERO
542:   100                CONTINUE
543:                      X( J ) = ONE
544:                      XJ = ONE
545:                      SCALE = ZERO
546:                      XMAX = ZERO
547:                   END IF
548:   105          CONTINUE
549: *
550: *              Scale x if necessary to avoid overflow when adding a
551: *              multiple of column j of A.
552: *
553:                IF( XJ.GT.ONE ) THEN
554:                   REC = ONE / XJ
555:                   IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
556: *
557: *                    Scale x by 1/(2*abs(x(j))).
558: *
559:                      REC = REC*HALF
560:                      CALL CSSCAL( N, REC, X, 1 )
561:                      SCALE = SCALE*REC
562:                   END IF
563:                ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
564: *
565: *                 Scale x by 1/2.
566: *
567:                   CALL CSSCAL( N, HALF, X, 1 )
568:                   SCALE = SCALE*HALF
569:                END IF
570: *
571:                IF( UPPER ) THEN
572:                   IF( J.GT.1 ) THEN
573: *
574: *                    Compute the update
575: *                       x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
576: *
577:                      CALL CAXPY( J-1, -X( J )*TSCAL, A( 1, J ), 1, X,
578:      \$                           1 )
579:                      I = ICAMAX( J-1, X, 1 )
580:                      XMAX = CABS1( X( I ) )
581:                   END IF
582:                ELSE
583:                   IF( J.LT.N ) THEN
584: *
585: *                    Compute the update
586: *                       x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
587: *
588:                      CALL CAXPY( N-J, -X( J )*TSCAL, A( J+1, J ), 1,
589:      \$                           X( J+1 ), 1 )
590:                      I = J + ICAMAX( N-J, X( J+1 ), 1 )
591:                      XMAX = CABS1( X( I ) )
592:                   END IF
593:                END IF
594:   110       CONTINUE
595: *
596:          ELSE IF( LSAME( TRANS, 'T' ) ) THEN
597: *
598: *           Solve A**T * x = b
599: *
600:             DO 150 J = JFIRST, JLAST, JINC
601: *
602: *              Compute x(j) = b(j) - sum A(k,j)*x(k).
603: *                                    k<>j
604: *
605:                XJ = CABS1( X( J ) )
606:                USCAL = TSCAL
607:                REC = ONE / MAX( XMAX, ONE )
608:                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
609: *
610: *                 If x(j) could overflow, scale x by 1/(2*XMAX).
611: *
612:                   REC = REC*HALF
613:                   IF( NOUNIT ) THEN
614:                      TJJS = A( J, J )*TSCAL
615:                   ELSE
616:                      TJJS = TSCAL
617:                   END IF
618:                      TJJ = CABS1( TJJS )
619:                      IF( TJJ.GT.ONE ) THEN
620: *
621: *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
622: *
623:                         REC = MIN( ONE, REC*TJJ )
624:                         USCAL = CLADIV( USCAL, TJJS )
625:                      END IF
626:                   IF( REC.LT.ONE ) THEN
627:                      CALL CSSCAL( N, REC, X, 1 )
628:                      SCALE = SCALE*REC
629:                      XMAX = XMAX*REC
630:                   END IF
631:                END IF
632: *
633:                CSUMJ = ZERO
634:                IF( USCAL.EQ.CMPLX( ONE ) ) THEN
635: *
636: *                 If the scaling needed for A in the dot product is 1,
637: *                 call CDOTU to perform the dot product.
638: *
639:                   IF( UPPER ) THEN
640:                      CSUMJ = CDOTU( J-1, A( 1, J ), 1, X, 1 )
641:                   ELSE IF( J.LT.N ) THEN
642:                      CSUMJ = CDOTU( N-J, A( J+1, J ), 1, X( J+1 ), 1 )
643:                   END IF
644:                ELSE
645: *
646: *                 Otherwise, use in-line code for the dot product.
647: *
648:                   IF( UPPER ) THEN
649:                      DO 120 I = 1, J - 1
650:                         CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
651:   120                CONTINUE
652:                   ELSE IF( J.LT.N ) THEN
653:                      DO 130 I = J + 1, N
654:                         CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
655:   130                CONTINUE
656:                   END IF
657:                END IF
658: *
659:                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
660: *
661: *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
662: *                 was not used to scale the dotproduct.
663: *
664:                   X( J ) = X( J ) - CSUMJ
665:                   XJ = CABS1( X( J ) )
666:                   IF( NOUNIT ) THEN
667:                      TJJS = A( J, J )*TSCAL
668:                   ELSE
669:                      TJJS = TSCAL
670:                      IF( TSCAL.EQ.ONE )
671:      \$                  GO TO 145
672:                   END IF
673: *
674: *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
675: *
676:                      TJJ = CABS1( TJJS )
677:                      IF( TJJ.GT.SMLNUM ) THEN
678: *
679: *                       abs(A(j,j)) > SMLNUM:
680: *
681:                         IF( TJJ.LT.ONE ) THEN
682:                            IF( XJ.GT.TJJ*BIGNUM ) THEN
683: *
684: *                             Scale X by 1/abs(x(j)).
685: *
686:                               REC = ONE / XJ
687:                               CALL CSSCAL( N, REC, X, 1 )
688:                               SCALE = SCALE*REC
689:                               XMAX = XMAX*REC
690:                            END IF
691:                         END IF
692:                         X( J ) = CLADIV( X( J ), TJJS )
693:                      ELSE IF( TJJ.GT.ZERO ) THEN
694: *
695: *                       0 < abs(A(j,j)) <= SMLNUM:
696: *
697:                         IF( XJ.GT.TJJ*BIGNUM ) THEN
698: *
699: *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
700: *
701:                            REC = ( TJJ*BIGNUM ) / XJ
702:                            CALL CSSCAL( N, REC, X, 1 )
703:                            SCALE = SCALE*REC
704:                            XMAX = XMAX*REC
705:                         END IF
706:                         X( J ) = CLADIV( X( J ), TJJS )
707:                      ELSE
708: *
709: *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
710: *                       scale = 0 and compute a solution to A**T *x = 0.
711: *
712:                         DO 140 I = 1, N
713:                            X( I ) = ZERO
714:   140                   CONTINUE
715:                         X( J ) = ONE
716:                         SCALE = ZERO
717:                         XMAX = ZERO
718:                      END IF
719:   145             CONTINUE
720:                ELSE
721: *
722: *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
723: *                 product has already been divided by 1/A(j,j).
724: *
725:                   X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
726:                END IF
727:                XMAX = MAX( XMAX, CABS1( X( J ) ) )
728:   150       CONTINUE
729: *
730:          ELSE
731: *
732: *           Solve A**H * x = b
733: *
734:             DO 190 J = JFIRST, JLAST, JINC
735: *
736: *              Compute x(j) = b(j) - sum A(k,j)*x(k).
737: *                                    k<>j
738: *
739:                XJ = CABS1( X( J ) )
740:                USCAL = TSCAL
741:                REC = ONE / MAX( XMAX, ONE )
742:                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
743: *
744: *                 If x(j) could overflow, scale x by 1/(2*XMAX).
745: *
746:                   REC = REC*HALF
747:                   IF( NOUNIT ) THEN
748:                      TJJS = CONJG( A( J, J ) )*TSCAL
749:                   ELSE
750:                      TJJS = TSCAL
751:                   END IF
752:                      TJJ = CABS1( TJJS )
753:                      IF( TJJ.GT.ONE ) THEN
754: *
755: *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
756: *
757:                         REC = MIN( ONE, REC*TJJ )
758:                         USCAL = CLADIV( USCAL, TJJS )
759:                      END IF
760:                   IF( REC.LT.ONE ) THEN
761:                      CALL CSSCAL( N, REC, X, 1 )
762:                      SCALE = SCALE*REC
763:                      XMAX = XMAX*REC
764:                   END IF
765:                END IF
766: *
767:                CSUMJ = ZERO
768:                IF( USCAL.EQ.CMPLX( ONE ) ) THEN
769: *
770: *                 If the scaling needed for A in the dot product is 1,
771: *                 call CDOTC to perform the dot product.
772: *
773:                   IF( UPPER ) THEN
774:                      CSUMJ = CDOTC( J-1, A( 1, J ), 1, X, 1 )
775:                   ELSE IF( J.LT.N ) THEN
776:                      CSUMJ = CDOTC( N-J, A( J+1, J ), 1, X( J+1 ), 1 )
777:                   END IF
778:                ELSE
779: *
780: *                 Otherwise, use in-line code for the dot product.
781: *
782:                   IF( UPPER ) THEN
783:                      DO 160 I = 1, J - 1
784:                         CSUMJ = CSUMJ + ( CONJG( A( I, J ) )*USCAL )*
785:      \$                          X( I )
786:   160                CONTINUE
787:                   ELSE IF( J.LT.N ) THEN
788:                      DO 170 I = J + 1, N
789:                         CSUMJ = CSUMJ + ( CONJG( A( I, J ) )*USCAL )*
790:      \$                          X( I )
791:   170                CONTINUE
792:                   END IF
793:                END IF
794: *
795:                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
796: *
797: *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
798: *                 was not used to scale the dotproduct.
799: *
800:                   X( J ) = X( J ) - CSUMJ
801:                   XJ = CABS1( X( J ) )
802:                   IF( NOUNIT ) THEN
803:                      TJJS = CONJG( A( J, J ) )*TSCAL
804:                   ELSE
805:                      TJJS = TSCAL
806:                      IF( TSCAL.EQ.ONE )
807:      \$                  GO TO 185
808:                   END IF
809: *
810: *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
811: *
812:                      TJJ = CABS1( TJJS )
813:                      IF( TJJ.GT.SMLNUM ) THEN
814: *
815: *                       abs(A(j,j)) > SMLNUM:
816: *
817:                         IF( TJJ.LT.ONE ) THEN
818:                            IF( XJ.GT.TJJ*BIGNUM ) THEN
819: *
820: *                             Scale X by 1/abs(x(j)).
821: *
822:                               REC = ONE / XJ
823:                               CALL CSSCAL( N, REC, X, 1 )
824:                               SCALE = SCALE*REC
825:                               XMAX = XMAX*REC
826:                            END IF
827:                         END IF
828:                         X( J ) = CLADIV( X( J ), TJJS )
829:                      ELSE IF( TJJ.GT.ZERO ) THEN
830: *
831: *                       0 < abs(A(j,j)) <= SMLNUM:
832: *
833:                         IF( XJ.GT.TJJ*BIGNUM ) THEN
834: *
835: *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
836: *
837:                            REC = ( TJJ*BIGNUM ) / XJ
838:                            CALL CSSCAL( N, REC, X, 1 )
839:                            SCALE = SCALE*REC
840:                            XMAX = XMAX*REC
841:                         END IF
842:                         X( J ) = CLADIV( X( J ), TJJS )
843:                      ELSE
844: *
845: *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
846: *                       scale = 0 and compute a solution to A**H *x = 0.
847: *
848:                         DO 180 I = 1, N
849:                            X( I ) = ZERO
850:   180                   CONTINUE
851:                         X( J ) = ONE
852:                         SCALE = ZERO
853:                         XMAX = ZERO
854:                      END IF
855:   185             CONTINUE
856:                ELSE
857: *
858: *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
859: *                 product has already been divided by 1/A(j,j).
860: *
861:                   X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
862:                END IF
863:                XMAX = MAX( XMAX, CABS1( X( J ) ) )
864:   190       CONTINUE
865:          END IF
866:          SCALE = SCALE / TSCAL
867:       END IF
868: *
869: *     Scale the column norms by 1/TSCAL for return.
870: *
871:       IF( TSCAL.NE.ONE ) THEN
872:          CALL SSCAL( N, ONE / TSCAL, CNORM, 1 )
873:       END IF
874: *
875:       RETURN
876: *
877: *     End of CLATRS
878: *
879:       END
880: ```