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