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: *> \date September 2012
155: *
156: *> \ingroup complex16OTHERauxiliary
157: *
158: *> \par Further Details:
159: * =====================
160: *>
161: *> \verbatim
162: *>
163: *> A rough bound on x is computed; if that is less than overflow, ZTPSV
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 ZTPSV 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 ZTPSV if 1/M(n) and 1/G(n) are both greater
227: *> than max(underflow, 1/overflow).
228: *> \endverbatim
229: *>
230: * =====================================================================
231: SUBROUTINE ZLATPS( 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: DOUBLE PRECISION SCALE
243: * ..
244: * .. Array Arguments ..
245: DOUBLE PRECISION CNORM( * )
246: COMPLEX*16 AP( * ), X( * )
247: * ..
248: *
249: * =====================================================================
250: *
251: * .. Parameters ..
252: DOUBLE PRECISION ZERO, HALF, ONE, TWO
253: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
254: $ TWO = 2.0D+0 )
255: * ..
256: * .. Local Scalars ..
257: LOGICAL NOTRAN, NOUNIT, UPPER
258: INTEGER I, IMAX, IP, J, JFIRST, JINC, JLAST, JLEN
259: DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
260: $ XBND, XJ, XMAX
261: COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
262: * ..
263: * .. External Functions ..
264: LOGICAL LSAME
265: INTEGER IDAMAX, IZAMAX
266: DOUBLE PRECISION DLAMCH, DZASUM
267: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
268: EXTERNAL LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC,
269: $ ZDOTU, ZLADIV
270: * ..
271: * .. External Subroutines ..
272: EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTPSV
273: * ..
274: * .. Intrinsic Functions ..
275: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
276: * ..
277: * .. Statement Functions ..
278: DOUBLE PRECISION CABS1, CABS2
279: * ..
280: * .. Statement Function definitions ..
281: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
282: CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) +
283: $ ABS( DIMAG( ZDUM ) / 2.D0 )
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( 'ZLATPS', -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 = DLAMCH( 'Safe minimum' )
320: BIGNUM = ONE / SMLNUM
321: CALL DLABAD( SMLNUM, BIGNUM )
322: SMLNUM = SMLNUM / DLAMCH( '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 ) = DZASUM( 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 ) = DZASUM( 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 = IDAMAX( 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 DSCAL( 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 ZTPSV 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 ZTPSV( 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 ZDSCAL( 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 120 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 110
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 ZDSCAL( N, REC, X, 1 )
595: SCALE = SCALE*REC
596: XMAX = XMAX*REC
597: END IF
598: END IF
599: X( J ) = ZLADIV( 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 ZDSCAL( N, REC, X, 1 )
619: SCALE = SCALE*REC
620: XMAX = XMAX*REC
621: END IF
622: X( J ) = ZLADIV( 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: 110 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 ZDSCAL( 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 ZDSCAL( 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 ZAXPY( J-1, -X( J )*TSCAL, AP( IP-J+1 ), 1, X,
667: $ 1 )
668: I = IZAMAX( 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 ZAXPY( N-J, -X( J )*TSCAL, AP( IP+1 ), 1,
679: $ X( J+1 ), 1 )
680: I = J + IZAMAX( N-J, X( J+1 ), 1 )
681: XMAX = CABS1( X( I ) )
682: END IF
683: IP = IP + N - J + 1
684: END IF
685: 120 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 170 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 = ZLADIV( USCAL, TJJS )
718: END IF
719: IF( REC.LT.ONE ) THEN
720: CALL ZDSCAL( 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.DCMPLX( ONE ) ) THEN
728: *
729: * If the scaling needed for A in the dot product is 1,
730: * call ZDOTU to perform the dot product.
731: *
732: IF( UPPER ) THEN
733: CSUMJ = ZDOTU( J-1, AP( IP-J+1 ), 1, X, 1 )
734: ELSE IF( J.LT.N ) THEN
735: CSUMJ = ZDOTU( 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 130 I = 1, J - 1
743: CSUMJ = CSUMJ + ( AP( IP-J+I )*USCAL )*X( I )
744: 130 CONTINUE
745: ELSE IF( J.LT.N ) THEN
746: DO 140 I = 1, N - J
747: CSUMJ = CSUMJ + ( AP( IP+I )*USCAL )*X( J+I )
748: 140 CONTINUE
749: END IF
750: END IF
751: *
752: IF( USCAL.EQ.DCMPLX( 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 160
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 ZDSCAL( N, REC, X, 1 )
781: SCALE = SCALE*REC
782: XMAX = XMAX*REC
783: END IF
784: END IF
785: X( J ) = ZLADIV( 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 ZDSCAL( N, REC, X, 1 )
796: SCALE = SCALE*REC
797: XMAX = XMAX*REC
798: END IF
799: X( J ) = ZLADIV( 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 150 I = 1, N
806: X( I ) = ZERO
807: 150 CONTINUE
808: X( J ) = ONE
809: SCALE = ZERO
810: XMAX = ZERO
811: END IF
812: 160 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 ) = ZLADIV( X( J ), TJJS ) - CSUMJ
819: END IF
820: XMAX = MAX( XMAX, CABS1( X( J ) ) )
821: JLEN = JLEN + 1
822: IP = IP + JINC*JLEN
823: 170 CONTINUE
824: *
825: ELSE
826: *
827: * Solve A**H * x = b
828: *
829: IP = JFIRST*( JFIRST+1 ) / 2
830: JLEN = 1
831: DO 220 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 = DCONJG( 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 = ZLADIV( USCAL, TJJS )
856: END IF
857: IF( REC.LT.ONE ) THEN
858: CALL ZDSCAL( 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.DCMPLX( ONE ) ) THEN
866: *
867: * If the scaling needed for A in the dot product is 1,
868: * call ZDOTC to perform the dot product.
869: *
870: IF( UPPER ) THEN
871: CSUMJ = ZDOTC( J-1, AP( IP-J+1 ), 1, X, 1 )
872: ELSE IF( J.LT.N ) THEN
873: CSUMJ = ZDOTC( 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 180 I = 1, J - 1
881: CSUMJ = CSUMJ + ( DCONJG( AP( IP-J+I ) )*USCAL )
882: $ *X( I )
883: 180 CONTINUE
884: ELSE IF( J.LT.N ) THEN
885: DO 190 I = 1, N - J
886: CSUMJ = CSUMJ + ( DCONJG( AP( IP+I ) )*USCAL )*
887: $ X( J+I )
888: 190 CONTINUE
889: END IF
890: END IF
891: *
892: IF( USCAL.EQ.DCMPLX( 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 = DCONJG( AP( IP ) )*TSCAL
904: ELSE
905: TJJS = TSCAL
906: IF( TSCAL.EQ.ONE )
907: $ GO TO 210
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 ZDSCAL( N, REC, X, 1 )
921: SCALE = SCALE*REC
922: XMAX = XMAX*REC
923: END IF
924: END IF
925: X( J ) = ZLADIV( 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 ZDSCAL( N, REC, X, 1 )
936: SCALE = SCALE*REC
937: XMAX = XMAX*REC
938: END IF
939: X( J ) = ZLADIV( 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 200 I = 1, N
946: X( I ) = ZERO
947: 200 CONTINUE
948: X( J ) = ONE
949: SCALE = ZERO
950: XMAX = ZERO
951: END IF
952: 210 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 ) = ZLADIV( X( J ), TJJS ) - CSUMJ
959: END IF
960: XMAX = MAX( XMAX, CABS1( X( J ) ) )
961: JLEN = JLEN + 1
962: IP = IP + JINC*JLEN
963: 220 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 DSCAL( N, ONE / TSCAL, CNORM, 1 )
972: END IF
973: *
974: RETURN
975: *
976: * End of ZLATPS
977: *
978: 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