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