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