1: SUBROUTINE ZLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
2: $ CNORM, INFO )
3: *
4: * -- LAPACK auxiliary routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER DIAG, NORMIN, TRANS, UPLO
11: INTEGER INFO, N
12: DOUBLE PRECISION SCALE
13: * ..
14: * .. Array Arguments ..
15: DOUBLE PRECISION CNORM( * )
16: COMPLEX*16 AP( * ), X( * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * ZLATPS solves one of the triangular systems
23: *
24: * A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
25: *
26: * with scaling to prevent overflow, where A is an upper or lower
27: * triangular matrix stored in packed form. Here A**T denotes the
28: * transpose of A, A**H denotes the conjugate transpose of A, x and b
29: * are n-element vectors, and s is a scaling factor, usually less than
30: * or equal to 1, chosen so that the components of x will be less than
31: * the overflow threshold. If the unscaled problem will not cause
32: * overflow, the Level 2 BLAS routine ZTPSV is called. If the matrix A
33: * is singular (A(j,j) = 0 for some j), then s is set to 0 and a
34: * non-trivial solution to A*x = 0 is returned.
35: *
36: * Arguments
37: * =========
38: *
39: * UPLO (input) CHARACTER*1
40: * Specifies whether the matrix A is upper or lower triangular.
41: * = 'U': Upper triangular
42: * = 'L': Lower triangular
43: *
44: * TRANS (input) CHARACTER*1
45: * Specifies the operation applied to A.
46: * = 'N': Solve A * x = s*b (No transpose)
47: * = 'T': Solve A**T * x = s*b (Transpose)
48: * = 'C': Solve A**H * x = s*b (Conjugate transpose)
49: *
50: * DIAG (input) CHARACTER*1
51: * Specifies whether or not the matrix A is unit triangular.
52: * = 'N': Non-unit triangular
53: * = 'U': Unit triangular
54: *
55: * NORMIN (input) CHARACTER*1
56: * Specifies whether CNORM has been set or not.
57: * = 'Y': CNORM contains the column norms on entry
58: * = 'N': CNORM is not set on entry. On exit, the norms will
59: * be computed and stored in CNORM.
60: *
61: * N (input) INTEGER
62: * The order of the matrix A. N >= 0.
63: *
64: * AP (input) COMPLEX*16 array, dimension (N*(N+1)/2)
65: * The upper or lower triangular matrix A, packed columnwise in
66: * a linear array. The j-th column of A is stored in the array
67: * AP as follows:
68: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
69: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
70: *
71: * X (input/output) COMPLEX*16 array, dimension (N)
72: * On entry, the right hand side b of the triangular system.
73: * On exit, X is overwritten by the solution vector x.
74: *
75: * SCALE (output) DOUBLE PRECISION
76: * The scaling factor s for the triangular system
77: * A * x = s*b, A**T * x = s*b, or A**H * x = s*b.
78: * If SCALE = 0, the matrix A is singular or badly scaled, and
79: * the vector x is an exact or approximate solution to A*x = 0.
80: *
81: * CNORM (input or output) DOUBLE PRECISION array, dimension (N)
82: *
83: * If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
84: * contains the norm of the off-diagonal part of the j-th column
85: * of A. If TRANS = 'N', CNORM(j) must be greater than or equal
86: * to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
87: * must be greater than or equal to the 1-norm.
88: *
89: * If NORMIN = 'N', CNORM is an output argument and CNORM(j)
90: * returns the 1-norm of the offdiagonal part of the j-th column
91: * of A.
92: *
93: * INFO (output) INTEGER
94: * = 0: successful exit
95: * < 0: if INFO = -k, the k-th argument had an illegal value
96: *
97: * Further Details
98: * ======= =======
99: *
100: * A rough bound on x is computed; if that is less than overflow, ZTPSV
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 ZTPSV 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 ZTPSV if 1/M(n) and 1/G(n) are both greater
164: * than max(underflow, 1/overflow).
165: *
166: * =====================================================================
167: *
168: * .. Parameters ..
169: DOUBLE PRECISION ZERO, HALF, ONE, TWO
170: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
171: $ TWO = 2.0D+0 )
172: * ..
173: * .. Local Scalars ..
174: LOGICAL NOTRAN, NOUNIT, UPPER
175: INTEGER I, IMAX, IP, J, JFIRST, JINC, JLAST, JLEN
176: DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
177: $ XBND, XJ, XMAX
178: COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
179: * ..
180: * .. External Functions ..
181: LOGICAL LSAME
182: INTEGER IDAMAX, IZAMAX
183: DOUBLE PRECISION DLAMCH, DZASUM
184: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
185: EXTERNAL LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC,
186: $ ZDOTU, ZLADIV
187: * ..
188: * .. External Subroutines ..
189: EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTPSV
190: * ..
191: * .. Intrinsic Functions ..
192: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
193: * ..
194: * .. Statement Functions ..
195: DOUBLE PRECISION CABS1, CABS2
196: * ..
197: * .. Statement Function definitions ..
198: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
199: CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) +
200: $ ABS( DIMAG( ZDUM ) / 2.D0 )
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( 'ZLATPS', -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 = DLAMCH( 'Safe minimum' )
237: BIGNUM = ONE / SMLNUM
238: CALL DLABAD( SMLNUM, BIGNUM )
239: SMLNUM = SMLNUM / DLAMCH( '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 ) = DZASUM( 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 ) = DZASUM( 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 = IDAMAX( 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 DSCAL( 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 ZTPSV 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 ZTPSV( 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 ZDSCAL( 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 120 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 110
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 ZDSCAL( N, REC, X, 1 )
512: SCALE = SCALE*REC
513: XMAX = XMAX*REC
514: END IF
515: END IF
516: X( J ) = ZLADIV( 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 ZDSCAL( N, REC, X, 1 )
536: SCALE = SCALE*REC
537: XMAX = XMAX*REC
538: END IF
539: X( J ) = ZLADIV( 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: 110 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 ZDSCAL( 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 ZDSCAL( 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 ZAXPY( J-1, -X( J )*TSCAL, AP( IP-J+1 ), 1, X,
584: $ 1 )
585: I = IZAMAX( 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 ZAXPY( N-J, -X( J )*TSCAL, AP( IP+1 ), 1,
596: $ X( J+1 ), 1 )
597: I = J + IZAMAX( N-J, X( J+1 ), 1 )
598: XMAX = CABS1( X( I ) )
599: END IF
600: IP = IP + N - J + 1
601: END IF
602: 120 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 170 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 = ZLADIV( USCAL, TJJS )
635: END IF
636: IF( REC.LT.ONE ) THEN
637: CALL ZDSCAL( 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.DCMPLX( ONE ) ) THEN
645: *
646: * If the scaling needed for A in the dot product is 1,
647: * call ZDOTU to perform the dot product.
648: *
649: IF( UPPER ) THEN
650: CSUMJ = ZDOTU( J-1, AP( IP-J+1 ), 1, X, 1 )
651: ELSE IF( J.LT.N ) THEN
652: CSUMJ = ZDOTU( 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 130 I = 1, J - 1
660: CSUMJ = CSUMJ + ( AP( IP-J+I )*USCAL )*X( I )
661: 130 CONTINUE
662: ELSE IF( J.LT.N ) THEN
663: DO 140 I = 1, N - J
664: CSUMJ = CSUMJ + ( AP( IP+I )*USCAL )*X( J+I )
665: 140 CONTINUE
666: END IF
667: END IF
668: *
669: IF( USCAL.EQ.DCMPLX( 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 160
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 ZDSCAL( N, REC, X, 1 )
698: SCALE = SCALE*REC
699: XMAX = XMAX*REC
700: END IF
701: END IF
702: X( J ) = ZLADIV( 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 ZDSCAL( N, REC, X, 1 )
713: SCALE = SCALE*REC
714: XMAX = XMAX*REC
715: END IF
716: X( J ) = ZLADIV( 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 150 I = 1, N
723: X( I ) = ZERO
724: 150 CONTINUE
725: X( J ) = ONE
726: SCALE = ZERO
727: XMAX = ZERO
728: END IF
729: 160 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 ) = ZLADIV( X( J ), TJJS ) - CSUMJ
736: END IF
737: XMAX = MAX( XMAX, CABS1( X( J ) ) )
738: JLEN = JLEN + 1
739: IP = IP + JINC*JLEN
740: 170 CONTINUE
741: *
742: ELSE
743: *
744: * Solve A**H * x = b
745: *
746: IP = JFIRST*( JFIRST+1 ) / 2
747: JLEN = 1
748: DO 220 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 = DCONJG( 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 = ZLADIV( USCAL, TJJS )
773: END IF
774: IF( REC.LT.ONE ) THEN
775: CALL ZDSCAL( 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.DCMPLX( ONE ) ) THEN
783: *
784: * If the scaling needed for A in the dot product is 1,
785: * call ZDOTC to perform the dot product.
786: *
787: IF( UPPER ) THEN
788: CSUMJ = ZDOTC( J-1, AP( IP-J+1 ), 1, X, 1 )
789: ELSE IF( J.LT.N ) THEN
790: CSUMJ = ZDOTC( 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 180 I = 1, J - 1
798: CSUMJ = CSUMJ + ( DCONJG( AP( IP-J+I ) )*USCAL )
799: $ *X( I )
800: 180 CONTINUE
801: ELSE IF( J.LT.N ) THEN
802: DO 190 I = 1, N - J
803: CSUMJ = CSUMJ + ( DCONJG( AP( IP+I ) )*USCAL )*
804: $ X( J+I )
805: 190 CONTINUE
806: END IF
807: END IF
808: *
809: IF( USCAL.EQ.DCMPLX( 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 = DCONJG( AP( IP ) )*TSCAL
821: ELSE
822: TJJS = TSCAL
823: IF( TSCAL.EQ.ONE )
824: $ GO TO 210
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 ZDSCAL( N, REC, X, 1 )
838: SCALE = SCALE*REC
839: XMAX = XMAX*REC
840: END IF
841: END IF
842: X( J ) = ZLADIV( 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 ZDSCAL( N, REC, X, 1 )
853: SCALE = SCALE*REC
854: XMAX = XMAX*REC
855: END IF
856: X( J ) = ZLADIV( 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 200 I = 1, N
863: X( I ) = ZERO
864: 200 CONTINUE
865: X( J ) = ONE
866: SCALE = ZERO
867: XMAX = ZERO
868: END IF
869: 210 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 ) = ZLADIV( X( J ), TJJS ) - CSUMJ
876: END IF
877: XMAX = MAX( XMAX, CABS1( X( J ) ) )
878: JLEN = JLEN + 1
879: IP = IP + JINC*JLEN
880: 220 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 DSCAL( N, ONE / TSCAL, CNORM, 1 )
889: END IF
890: *
891: RETURN
892: *
893: * End of ZLATPS
894: *
895: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlatps.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack