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