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