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