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: