ENH: Improving the travis dashboard name
[platform/upstream/lapack.git] / SRC / clatps.f
1 *> \brief \b CLATPS solves a triangular system of equations with the matrix held in packed storage.
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CLATPS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatps.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatps.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatps.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE CLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
22 *                          CNORM, INFO )
23 *
24 *       .. Scalar Arguments ..
25 *       CHARACTER          DIAG, NORMIN, TRANS, UPLO
26 *       INTEGER            INFO, N
27 *       REAL               SCALE
28 *       ..
29 *       .. Array Arguments ..
30 *       REAL               CNORM( * )
31 *       COMPLEX            AP( * ), X( * )
32 *       ..
33 *
34 *
35 *> \par Purpose:
36 *  =============
37 *>
38 *> \verbatim
39 *>
40 *> CLATPS solves one of the triangular systems
41 *>
42 *>    A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,
43 *>
44 *> with scaling to prevent overflow, where A is an upper or lower
45 *> triangular matrix stored in packed form.  Here A**T denotes the
46 *> transpose of A, A**H denotes the conjugate transpose of A, x and b
47 *> are n-element vectors, and s is a scaling factor, usually less than
48 *> or equal to 1, chosen so that the components of x will be less than
49 *> the overflow threshold.  If the unscaled problem will not cause
50 *> overflow, the Level 2 BLAS routine CTPSV is called. If the matrix A
51 *> is singular (A(j,j) = 0 for some j), then s is set to 0 and a
52 *> non-trivial solution to A*x = 0 is returned.
53 *> \endverbatim
54 *
55 *  Arguments:
56 *  ==========
57 *
58 *> \param[in] UPLO
59 *> \verbatim
60 *>          UPLO is CHARACTER*1
61 *>          Specifies whether the matrix A is upper or lower triangular.
62 *>          = 'U':  Upper triangular
63 *>          = 'L':  Lower triangular
64 *> \endverbatim
65 *>
66 *> \param[in] TRANS
67 *> \verbatim
68 *>          TRANS is CHARACTER*1
69 *>          Specifies the operation applied to A.
70 *>          = 'N':  Solve A * x = s*b     (No transpose)
71 *>          = 'T':  Solve A**T * x = s*b  (Transpose)
72 *>          = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
73 *> \endverbatim
74 *>
75 *> \param[in] DIAG
76 *> \verbatim
77 *>          DIAG is CHARACTER*1
78 *>          Specifies whether or not the matrix A is unit triangular.
79 *>          = 'N':  Non-unit triangular
80 *>          = 'U':  Unit triangular
81 *> \endverbatim
82 *>
83 *> \param[in] NORMIN
84 *> \verbatim
85 *>          NORMIN is CHARACTER*1
86 *>          Specifies whether CNORM has been set or not.
87 *>          = 'Y':  CNORM contains the column norms on entry
88 *>          = 'N':  CNORM is not set on entry.  On exit, the norms will
89 *>                  be computed and stored in CNORM.
90 *> \endverbatim
91 *>
92 *> \param[in] N
93 *> \verbatim
94 *>          N is INTEGER
95 *>          The order of the matrix A.  N >= 0.
96 *> \endverbatim
97 *>
98 *> \param[in] AP
99 *> \verbatim
100 *>          AP is COMPLEX array, dimension (N*(N+1)/2)
101 *>          The upper or lower triangular matrix A, packed columnwise in
102 *>          a linear array.  The j-th column of A is stored in the array
103 *>          AP as follows:
104 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
105 *>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
106 *> \endverbatim
107 *>
108 *> \param[in,out] X
109 *> \verbatim
110 *>          X is COMPLEX array, dimension (N)
111 *>          On entry, the right hand side b of the triangular system.
112 *>          On exit, X is overwritten by the solution vector x.
113 *> \endverbatim
114 *>
115 *> \param[out] SCALE
116 *> \verbatim
117 *>          SCALE is REAL
118 *>          The scaling factor s for the triangular system
119 *>             A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b.
120 *>          If SCALE = 0, the matrix A is singular or badly scaled, and
121 *>          the vector x is an exact or approximate solution to A*x = 0.
122 *> \endverbatim
123 *>
124 *> \param[in,out] CNORM
125 *> \verbatim
126 *>          CNORM is REAL array, dimension (N)
127 *>
128 *>          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
129 *>          contains the norm of the off-diagonal part of the j-th column
130 *>          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
131 *>          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
132 *>          must be greater than or equal to the 1-norm.
133 *>
134 *>          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
135 *>          returns the 1-norm of the offdiagonal part of the j-th column
136 *>          of A.
137 *> \endverbatim
138 *>
139 *> \param[out] INFO
140 *> \verbatim
141 *>          INFO is INTEGER
142 *>          = 0:  successful exit
143 *>          < 0:  if INFO = -k, the k-th argument had an illegal value
144 *> \endverbatim
145 *
146 *  Authors:
147 *  ========
148 *
149 *> \author Univ. of Tennessee
150 *> \author Univ. of California Berkeley
151 *> \author Univ. of Colorado Denver
152 *> \author NAG Ltd.
153 *
154 *> \date September 2012
155 *
156 *> \ingroup complexOTHERauxiliary
157 *
158 *> \par Further Details:
159 *  =====================
160 *>
161 *> \verbatim
162 *>
163 *>  A rough bound on x is computed; if that is less than overflow, CTPSV
164 *>  is called, otherwise, specific code is used which checks for possible
165 *>  overflow or divide-by-zero at every operation.
166 *>
167 *>  A columnwise scheme is used for solving A*x = b.  The basic algorithm
168 *>  if A is lower triangular is
169 *>
170 *>       x[1:n] := b[1:n]
171 *>       for j = 1, ..., n
172 *>            x(j) := x(j) / A(j,j)
173 *>            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
174 *>       end
175 *>
176 *>  Define bounds on the components of x after j iterations of the loop:
177 *>     M(j) = bound on x[1:j]
178 *>     G(j) = bound on x[j+1:n]
179 *>  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
180 *>
181 *>  Then for iteration j+1 we have
182 *>     M(j+1) <= G(j) / | A(j+1,j+1) |
183 *>     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
184 *>            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
185 *>
186 *>  where CNORM(j+1) is greater than or equal to the infinity-norm of
187 *>  column j+1 of A, not counting the diagonal.  Hence
188 *>
189 *>     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
190 *>                  1<=i<=j
191 *>  and
192 *>
193 *>     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
194 *>                                   1<=i< j
195 *>
196 *>  Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTPSV if the
197 *>  reciprocal of the largest M(j), j=1,..,n, is larger than
198 *>  max(underflow, 1/overflow).
199 *>
200 *>  The bound on x(j) is also used to determine when a step in the
201 *>  columnwise method can be performed without fear of overflow.  If
202 *>  the computed bound is greater than a large constant, x is scaled to
203 *>  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
204 *>  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
205 *>
206 *>  Similarly, a row-wise scheme is used to solve A**T *x = b  or
207 *>  A**H *x = b.  The basic algorithm for A upper triangular is
208 *>
209 *>       for j = 1, ..., n
210 *>            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
211 *>       end
212 *>
213 *>  We simultaneously compute two bounds
214 *>       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
215 *>       M(j) = bound on x(i), 1<=i<=j
216 *>
217 *>  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
218 *>  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
219 *>  Then the bound on x(j) is
220 *>
221 *>       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
222 *>
223 *>            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
224 *>                      1<=i<=j
225 *>
226 *>  and we can safely call CTPSV if 1/M(n) and 1/G(n) are both greater
227 *>  than max(underflow, 1/overflow).
228 *> \endverbatim
229 *>
230 *  =====================================================================
231       SUBROUTINE CLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
232      $                   CNORM, INFO )
233 *
234 *  -- LAPACK auxiliary routine (version 3.4.2) --
235 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
236 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
237 *     September 2012
238 *
239 *     .. Scalar Arguments ..
240       CHARACTER          DIAG, NORMIN, TRANS, UPLO
241       INTEGER            INFO, N
242       REAL               SCALE
243 *     ..
244 *     .. Array Arguments ..
245       REAL               CNORM( * )
246       COMPLEX            AP( * ), X( * )
247 *     ..
248 *
249 *  =====================================================================
250 *
251 *     .. Parameters ..
252       REAL               ZERO, HALF, ONE, TWO
253       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0,
254      $                   TWO = 2.0E+0 )
255 *     ..
256 *     .. Local Scalars ..
257       LOGICAL            NOTRAN, NOUNIT, UPPER
258       INTEGER            I, IMAX, IP, J, JFIRST, JINC, JLAST, JLEN
259       REAL               BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
260      $                   XBND, XJ, XMAX
261       COMPLEX            CSUMJ, TJJS, USCAL, ZDUM
262 *     ..
263 *     .. External Functions ..
264       LOGICAL            LSAME
265       INTEGER            ICAMAX, ISAMAX
266       REAL               SCASUM, SLAMCH
267       COMPLEX            CDOTC, CDOTU, CLADIV
268       EXTERNAL           LSAME, ICAMAX, ISAMAX, SCASUM, SLAMCH, CDOTC,
269      $                   CDOTU, CLADIV
270 *     ..
271 *     .. External Subroutines ..
272       EXTERNAL           CAXPY, CSSCAL, CTPSV, SLABAD, SSCAL, XERBLA
273 *     ..
274 *     .. Intrinsic Functions ..
275       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
276 *     ..
277 *     .. Statement Functions ..
278       REAL               CABS1, CABS2
279 *     ..
280 *     .. Statement Function definitions ..
281       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
282       CABS2( ZDUM ) = ABS( REAL( ZDUM ) / 2. ) +
283      $                ABS( AIMAG( ZDUM ) / 2. )
284 *     ..
285 *     .. Executable Statements ..
286 *
287       INFO = 0
288       UPPER = LSAME( UPLO, 'U' )
289       NOTRAN = LSAME( TRANS, 'N' )
290       NOUNIT = LSAME( DIAG, 'N' )
291 *
292 *     Test the input parameters.
293 *
294       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
295          INFO = -1
296       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
297      $         LSAME( TRANS, 'C' ) ) THEN
298          INFO = -2
299       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
300          INFO = -3
301       ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
302      $         LSAME( NORMIN, 'N' ) ) THEN
303          INFO = -4
304       ELSE IF( N.LT.0 ) THEN
305          INFO = -5
306       END IF
307       IF( INFO.NE.0 ) THEN
308          CALL XERBLA( 'CLATPS', -INFO )
309          RETURN
310       END IF
311 *
312 *     Quick return if possible
313 *
314       IF( N.EQ.0 )
315      $   RETURN
316 *
317 *     Determine machine dependent parameters to control overflow.
318 *
319       SMLNUM = SLAMCH( 'Safe minimum' )
320       BIGNUM = ONE / SMLNUM
321       CALL SLABAD( SMLNUM, BIGNUM )
322       SMLNUM = SMLNUM / SLAMCH( 'Precision' )
323       BIGNUM = ONE / SMLNUM
324       SCALE = ONE
325 *
326       IF( LSAME( NORMIN, 'N' ) ) THEN
327 *
328 *        Compute the 1-norm of each column, not including the diagonal.
329 *
330          IF( UPPER ) THEN
331 *
332 *           A is upper triangular.
333 *
334             IP = 1
335             DO 10 J = 1, N
336                CNORM( J ) = SCASUM( J-1, AP( IP ), 1 )
337                IP = IP + J
338    10       CONTINUE
339          ELSE
340 *
341 *           A is lower triangular.
342 *
343             IP = 1
344             DO 20 J = 1, N - 1
345                CNORM( J ) = SCASUM( N-J, AP( IP+1 ), 1 )
346                IP = IP + N - J + 1
347    20       CONTINUE
348             CNORM( N ) = ZERO
349          END IF
350       END IF
351 *
352 *     Scale the column norms by TSCAL if the maximum element in CNORM is
353 *     greater than BIGNUM/2.
354 *
355       IMAX = ISAMAX( N, CNORM, 1 )
356       TMAX = CNORM( IMAX )
357       IF( TMAX.LE.BIGNUM*HALF ) THEN
358          TSCAL = ONE
359       ELSE
360          TSCAL = HALF / ( SMLNUM*TMAX )
361          CALL SSCAL( N, TSCAL, CNORM, 1 )
362       END IF
363 *
364 *     Compute a bound on the computed solution vector to see if the
365 *     Level 2 BLAS routine CTPSV can be used.
366 *
367       XMAX = ZERO
368       DO 30 J = 1, N
369          XMAX = MAX( XMAX, CABS2( X( J ) ) )
370    30 CONTINUE
371       XBND = XMAX
372       IF( NOTRAN ) THEN
373 *
374 *        Compute the growth in A * x = b.
375 *
376          IF( UPPER ) THEN
377             JFIRST = N
378             JLAST = 1
379             JINC = -1
380          ELSE
381             JFIRST = 1
382             JLAST = N
383             JINC = 1
384          END IF
385 *
386          IF( TSCAL.NE.ONE ) THEN
387             GROW = ZERO
388             GO TO 60
389          END IF
390 *
391          IF( NOUNIT ) THEN
392 *
393 *           A is non-unit triangular.
394 *
395 *           Compute GROW = 1/G(j) and XBND = 1/M(j).
396 *           Initially, G(0) = max{x(i), i=1,...,n}.
397 *
398             GROW = HALF / MAX( XBND, SMLNUM )
399             XBND = GROW
400             IP = JFIRST*( JFIRST+1 ) / 2
401             JLEN = N
402             DO 40 J = JFIRST, JLAST, JINC
403 *
404 *              Exit the loop if the growth factor is too small.
405 *
406                IF( GROW.LE.SMLNUM )
407      $            GO TO 60
408 *
409                TJJS = AP( IP )
410                TJJ = CABS1( TJJS )
411 *
412                IF( TJJ.GE.SMLNUM ) THEN
413 *
414 *                 M(j) = G(j-1) / abs(A(j,j))
415 *
416                   XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
417                ELSE
418 *
419 *                 M(j) could overflow, set XBND to 0.
420 *
421                   XBND = ZERO
422                END IF
423 *
424                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
425 *
426 *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
427 *
428                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
429                ELSE
430 *
431 *                 G(j) could overflow, set GROW to 0.
432 *
433                   GROW = ZERO
434                END IF
435                IP = IP + JINC*JLEN
436                JLEN = JLEN - 1
437    40       CONTINUE
438             GROW = XBND
439          ELSE
440 *
441 *           A is unit triangular.
442 *
443 *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
444 *
445             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
446             DO 50 J = JFIRST, JLAST, JINC
447 *
448 *              Exit the loop if the growth factor is too small.
449 *
450                IF( GROW.LE.SMLNUM )
451      $            GO TO 60
452 *
453 *              G(j) = G(j-1)*( 1 + CNORM(j) )
454 *
455                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
456    50       CONTINUE
457          END IF
458    60    CONTINUE
459 *
460       ELSE
461 *
462 *        Compute the growth in A**T * x = b  or  A**H * x = b.
463 *
464          IF( UPPER ) THEN
465             JFIRST = 1
466             JLAST = N
467             JINC = 1
468          ELSE
469             JFIRST = N
470             JLAST = 1
471             JINC = -1
472          END IF
473 *
474          IF( TSCAL.NE.ONE ) THEN
475             GROW = ZERO
476             GO TO 90
477          END IF
478 *
479          IF( NOUNIT ) THEN
480 *
481 *           A is non-unit triangular.
482 *
483 *           Compute GROW = 1/G(j) and XBND = 1/M(j).
484 *           Initially, M(0) = max{x(i), i=1,...,n}.
485 *
486             GROW = HALF / MAX( XBND, SMLNUM )
487             XBND = GROW
488             IP = JFIRST*( JFIRST+1 ) / 2
489             JLEN = 1
490             DO 70 J = JFIRST, JLAST, JINC
491 *
492 *              Exit the loop if the growth factor is too small.
493 *
494                IF( GROW.LE.SMLNUM )
495      $            GO TO 90
496 *
497 *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
498 *
499                XJ = ONE + CNORM( J )
500                GROW = MIN( GROW, XBND / XJ )
501 *
502                TJJS = AP( IP )
503                TJJ = CABS1( TJJS )
504 *
505                IF( TJJ.GE.SMLNUM ) THEN
506 *
507 *                 M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
508 *
509                   IF( XJ.GT.TJJ )
510      $               XBND = XBND*( TJJ / XJ )
511                ELSE
512 *
513 *                 M(j) could overflow, set XBND to 0.
514 *
515                   XBND = ZERO
516                END IF
517                JLEN = JLEN + 1
518                IP = IP + JINC*JLEN
519    70       CONTINUE
520             GROW = MIN( GROW, XBND )
521          ELSE
522 *
523 *           A is unit triangular.
524 *
525 *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
526 *
527             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
528             DO 80 J = JFIRST, JLAST, JINC
529 *
530 *              Exit the loop if the growth factor is too small.
531 *
532                IF( GROW.LE.SMLNUM )
533      $            GO TO 90
534 *
535 *              G(j) = ( 1 + CNORM(j) )*G(j-1)
536 *
537                XJ = ONE + CNORM( J )
538                GROW = GROW / XJ
539    80       CONTINUE
540          END IF
541    90    CONTINUE
542       END IF
543 *
544       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
545 *
546 *        Use the Level 2 BLAS solve if the reciprocal of the bound on
547 *        elements of X is not too small.
548 *
549          CALL CTPSV( UPLO, TRANS, DIAG, N, AP, X, 1 )
550       ELSE
551 *
552 *        Use a Level 1 BLAS solve, scaling intermediate results.
553 *
554          IF( XMAX.GT.BIGNUM*HALF ) THEN
555 *
556 *           Scale X so that its components are less than or equal to
557 *           BIGNUM in absolute value.
558 *
559             SCALE = ( BIGNUM*HALF ) / XMAX
560             CALL CSSCAL( N, SCALE, X, 1 )
561             XMAX = BIGNUM
562          ELSE
563             XMAX = XMAX*TWO
564          END IF
565 *
566          IF( NOTRAN ) THEN
567 *
568 *           Solve A * x = b
569 *
570             IP = JFIRST*( JFIRST+1 ) / 2
571             DO 110 J = JFIRST, JLAST, JINC
572 *
573 *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
574 *
575                XJ = CABS1( X( J ) )
576                IF( NOUNIT ) THEN
577                   TJJS = AP( IP )*TSCAL
578                ELSE
579                   TJJS = TSCAL
580                   IF( TSCAL.EQ.ONE )
581      $               GO TO 105
582                END IF
583                   TJJ = CABS1( TJJS )
584                   IF( TJJ.GT.SMLNUM ) THEN
585 *
586 *                    abs(A(j,j)) > SMLNUM:
587 *
588                      IF( TJJ.LT.ONE ) THEN
589                         IF( XJ.GT.TJJ*BIGNUM ) THEN
590 *
591 *                          Scale x by 1/b(j).
592 *
593                            REC = ONE / XJ
594                            CALL CSSCAL( N, REC, X, 1 )
595                            SCALE = SCALE*REC
596                            XMAX = XMAX*REC
597                         END IF
598                      END IF
599                      X( J ) = CLADIV( X( J ), TJJS )
600                      XJ = CABS1( X( J ) )
601                   ELSE IF( TJJ.GT.ZERO ) THEN
602 *
603 *                    0 < abs(A(j,j)) <= SMLNUM:
604 *
605                      IF( XJ.GT.TJJ*BIGNUM ) THEN
606 *
607 *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
608 *                       to avoid overflow when dividing by A(j,j).
609 *
610                         REC = ( TJJ*BIGNUM ) / XJ
611                         IF( CNORM( J ).GT.ONE ) THEN
612 *
613 *                          Scale by 1/CNORM(j) to avoid overflow when
614 *                          multiplying x(j) times column j.
615 *
616                            REC = REC / CNORM( J )
617                         END IF
618                         CALL CSSCAL( N, REC, X, 1 )
619                         SCALE = SCALE*REC
620                         XMAX = XMAX*REC
621                      END IF
622                      X( J ) = CLADIV( X( J ), TJJS )
623                      XJ = CABS1( X( J ) )
624                   ELSE
625 *
626 *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
627 *                    scale = 0, and compute a solution to A*x = 0.
628 *
629                      DO 100 I = 1, N
630                         X( I ) = ZERO
631   100                CONTINUE
632                      X( J ) = ONE
633                      XJ = ONE
634                      SCALE = ZERO
635                      XMAX = ZERO
636                   END IF
637   105          CONTINUE
638 *
639 *              Scale x if necessary to avoid overflow when adding a
640 *              multiple of column j of A.
641 *
642                IF( XJ.GT.ONE ) THEN
643                   REC = ONE / XJ
644                   IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
645 *
646 *                    Scale x by 1/(2*abs(x(j))).
647 *
648                      REC = REC*HALF
649                      CALL CSSCAL( N, REC, X, 1 )
650                      SCALE = SCALE*REC
651                   END IF
652                ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
653 *
654 *                 Scale x by 1/2.
655 *
656                   CALL CSSCAL( N, HALF, X, 1 )
657                   SCALE = SCALE*HALF
658                END IF
659 *
660                IF( UPPER ) THEN
661                   IF( J.GT.1 ) THEN
662 *
663 *                    Compute the update
664 *                       x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
665 *
666                      CALL CAXPY( J-1, -X( J )*TSCAL, AP( IP-J+1 ), 1, X,
667      $                           1 )
668                      I = ICAMAX( J-1, X, 1 )
669                      XMAX = CABS1( X( I ) )
670                   END IF
671                   IP = IP - J
672                ELSE
673                   IF( J.LT.N ) THEN
674 *
675 *                    Compute the update
676 *                       x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
677 *
678                      CALL CAXPY( N-J, -X( J )*TSCAL, AP( IP+1 ), 1,
679      $                           X( J+1 ), 1 )
680                      I = J + ICAMAX( N-J, X( J+1 ), 1 )
681                      XMAX = CABS1( X( I ) )
682                   END IF
683                   IP = IP + N - J + 1
684                END IF
685   110       CONTINUE
686 *
687          ELSE IF( LSAME( TRANS, 'T' ) ) THEN
688 *
689 *           Solve A**T * x = b
690 *
691             IP = JFIRST*( JFIRST+1 ) / 2
692             JLEN = 1
693             DO 150 J = JFIRST, JLAST, JINC
694 *
695 *              Compute x(j) = b(j) - sum A(k,j)*x(k).
696 *                                    k<>j
697 *
698                XJ = CABS1( X( J ) )
699                USCAL = TSCAL
700                REC = ONE / MAX( XMAX, ONE )
701                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
702 *
703 *                 If x(j) could overflow, scale x by 1/(2*XMAX).
704 *
705                   REC = REC*HALF
706                   IF( NOUNIT ) THEN
707                      TJJS = AP( IP )*TSCAL
708                   ELSE
709                      TJJS = TSCAL
710                   END IF
711                      TJJ = CABS1( TJJS )
712                      IF( TJJ.GT.ONE ) THEN
713 *
714 *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
715 *
716                         REC = MIN( ONE, REC*TJJ )
717                         USCAL = CLADIV( USCAL, TJJS )
718                      END IF
719                   IF( REC.LT.ONE ) THEN
720                      CALL CSSCAL( N, REC, X, 1 )
721                      SCALE = SCALE*REC
722                      XMAX = XMAX*REC
723                   END IF
724                END IF
725 *
726                CSUMJ = ZERO
727                IF( USCAL.EQ.CMPLX( ONE ) ) THEN
728 *
729 *                 If the scaling needed for A in the dot product is 1,
730 *                 call CDOTU to perform the dot product.
731 *
732                   IF( UPPER ) THEN
733                      CSUMJ = CDOTU( J-1, AP( IP-J+1 ), 1, X, 1 )
734                   ELSE IF( J.LT.N ) THEN
735                      CSUMJ = CDOTU( N-J, AP( IP+1 ), 1, X( J+1 ), 1 )
736                   END IF
737                ELSE
738 *
739 *                 Otherwise, use in-line code for the dot product.
740 *
741                   IF( UPPER ) THEN
742                      DO 120 I = 1, J - 1
743                         CSUMJ = CSUMJ + ( AP( IP-J+I )*USCAL )*X( I )
744   120                CONTINUE
745                   ELSE IF( J.LT.N ) THEN
746                      DO 130 I = 1, N - J
747                         CSUMJ = CSUMJ + ( AP( IP+I )*USCAL )*X( J+I )
748   130                CONTINUE
749                   END IF
750                END IF
751 *
752                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
753 *
754 *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
755 *                 was not used to scale the dotproduct.
756 *
757                   X( J ) = X( J ) - CSUMJ
758                   XJ = CABS1( X( J ) )
759                   IF( NOUNIT ) THEN
760 *
761 *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
762 *
763                      TJJS = AP( IP )*TSCAL
764                   ELSE
765                      TJJS = TSCAL
766                      IF( TSCAL.EQ.ONE )
767      $                  GO TO 145
768                   END IF
769                      TJJ = CABS1( TJJS )
770                      IF( TJJ.GT.SMLNUM ) THEN
771 *
772 *                       abs(A(j,j)) > SMLNUM:
773 *
774                         IF( TJJ.LT.ONE ) THEN
775                            IF( XJ.GT.TJJ*BIGNUM ) THEN
776 *
777 *                             Scale X by 1/abs(x(j)).
778 *
779                               REC = ONE / XJ
780                               CALL CSSCAL( N, REC, X, 1 )
781                               SCALE = SCALE*REC
782                               XMAX = XMAX*REC
783                            END IF
784                         END IF
785                         X( J ) = CLADIV( X( J ), TJJS )
786                      ELSE IF( TJJ.GT.ZERO ) THEN
787 *
788 *                       0 < abs(A(j,j)) <= SMLNUM:
789 *
790                         IF( XJ.GT.TJJ*BIGNUM ) THEN
791 *
792 *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
793 *
794                            REC = ( TJJ*BIGNUM ) / XJ
795                            CALL CSSCAL( N, REC, X, 1 )
796                            SCALE = SCALE*REC
797                            XMAX = XMAX*REC
798                         END IF
799                         X( J ) = CLADIV( X( J ), TJJS )
800                      ELSE
801 *
802 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
803 *                       scale = 0 and compute a solution to A**T *x = 0.
804 *
805                         DO 140 I = 1, N
806                            X( I ) = ZERO
807   140                   CONTINUE
808                         X( J ) = ONE
809                         SCALE = ZERO
810                         XMAX = ZERO
811                      END IF
812   145             CONTINUE
813                ELSE
814 *
815 *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
816 *                 product has already been divided by 1/A(j,j).
817 *
818                   X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
819                END IF
820                XMAX = MAX( XMAX, CABS1( X( J ) ) )
821                JLEN = JLEN + 1
822                IP = IP + JINC*JLEN
823   150       CONTINUE
824 *
825          ELSE
826 *
827 *           Solve A**H * x = b
828 *
829             IP = JFIRST*( JFIRST+1 ) / 2
830             JLEN = 1
831             DO 190 J = JFIRST, JLAST, JINC
832 *
833 *              Compute x(j) = b(j) - sum A(k,j)*x(k).
834 *                                    k<>j
835 *
836                XJ = CABS1( X( J ) )
837                USCAL = TSCAL
838                REC = ONE / MAX( XMAX, ONE )
839                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
840 *
841 *                 If x(j) could overflow, scale x by 1/(2*XMAX).
842 *
843                   REC = REC*HALF
844                   IF( NOUNIT ) THEN
845                      TJJS = CONJG( AP( IP ) )*TSCAL
846                   ELSE
847                      TJJS = TSCAL
848                   END IF
849                      TJJ = CABS1( TJJS )
850                      IF( TJJ.GT.ONE ) THEN
851 *
852 *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
853 *
854                         REC = MIN( ONE, REC*TJJ )
855                         USCAL = CLADIV( USCAL, TJJS )
856                      END IF
857                   IF( REC.LT.ONE ) THEN
858                      CALL CSSCAL( N, REC, X, 1 )
859                      SCALE = SCALE*REC
860                      XMAX = XMAX*REC
861                   END IF
862                END IF
863 *
864                CSUMJ = ZERO
865                IF( USCAL.EQ.CMPLX( ONE ) ) THEN
866 *
867 *                 If the scaling needed for A in the dot product is 1,
868 *                 call CDOTC to perform the dot product.
869 *
870                   IF( UPPER ) THEN
871                      CSUMJ = CDOTC( J-1, AP( IP-J+1 ), 1, X, 1 )
872                   ELSE IF( J.LT.N ) THEN
873                      CSUMJ = CDOTC( N-J, AP( IP+1 ), 1, X( J+1 ), 1 )
874                   END IF
875                ELSE
876 *
877 *                 Otherwise, use in-line code for the dot product.
878 *
879                   IF( UPPER ) THEN
880                      DO 160 I = 1, J - 1
881                         CSUMJ = CSUMJ + ( CONJG( AP( IP-J+I ) )*USCAL )*
882      $                          X( I )
883   160                CONTINUE
884                   ELSE IF( J.LT.N ) THEN
885                      DO 170 I = 1, N - J
886                         CSUMJ = CSUMJ + ( CONJG( AP( IP+I ) )*USCAL )*
887      $                          X( J+I )
888   170                CONTINUE
889                   END IF
890                END IF
891 *
892                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
893 *
894 *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
895 *                 was not used to scale the dotproduct.
896 *
897                   X( J ) = X( J ) - CSUMJ
898                   XJ = CABS1( X( J ) )
899                   IF( NOUNIT ) THEN
900 *
901 *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
902 *
903                      TJJS = CONJG( AP( IP ) )*TSCAL
904                   ELSE
905                      TJJS = TSCAL
906                      IF( TSCAL.EQ.ONE )
907      $                  GO TO 185
908                   END IF
909                      TJJ = CABS1( TJJS )
910                      IF( TJJ.GT.SMLNUM ) THEN
911 *
912 *                       abs(A(j,j)) > SMLNUM:
913 *
914                         IF( TJJ.LT.ONE ) THEN
915                            IF( XJ.GT.TJJ*BIGNUM ) THEN
916 *
917 *                             Scale X by 1/abs(x(j)).
918 *
919                               REC = ONE / XJ
920                               CALL CSSCAL( N, REC, X, 1 )
921                               SCALE = SCALE*REC
922                               XMAX = XMAX*REC
923                            END IF
924                         END IF
925                         X( J ) = CLADIV( X( J ), TJJS )
926                      ELSE IF( TJJ.GT.ZERO ) THEN
927 *
928 *                       0 < abs(A(j,j)) <= SMLNUM:
929 *
930                         IF( XJ.GT.TJJ*BIGNUM ) THEN
931 *
932 *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
933 *
934                            REC = ( TJJ*BIGNUM ) / XJ
935                            CALL CSSCAL( N, REC, X, 1 )
936                            SCALE = SCALE*REC
937                            XMAX = XMAX*REC
938                         END IF
939                         X( J ) = CLADIV( X( J ), TJJS )
940                      ELSE
941 *
942 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
943 *                       scale = 0 and compute a solution to A**H *x = 0.
944 *
945                         DO 180 I = 1, N
946                            X( I ) = ZERO
947   180                   CONTINUE
948                         X( J ) = ONE
949                         SCALE = ZERO
950                         XMAX = ZERO
951                      END IF
952   185             CONTINUE
953                ELSE
954 *
955 *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
956 *                 product has already been divided by 1/A(j,j).
957 *
958                   X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
959                END IF
960                XMAX = MAX( XMAX, CABS1( X( J ) ) )
961                JLEN = JLEN + 1
962                IP = IP + JINC*JLEN
963   190       CONTINUE
964          END IF
965          SCALE = SCALE / TSCAL
966       END IF
967 *
968 *     Scale the column norms by 1/TSCAL for return.
969 *
970       IF( TSCAL.NE.ONE ) THEN
971          CALL SSCAL( N, ONE / TSCAL, CNORM, 1 )
972       END IF
973 *
974       RETURN
975 *
976 *     End of CLATPS
977 *
978       END