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**T*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**T 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**T* x = s*b (Transpose)
46 * = 'C': Solve A**T* 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**T* 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**T*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]**T * 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]**T * 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**T * 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**T * 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**T*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
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**T*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**T 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**T* x = s*b (Transpose)
46 * = 'C': Solve A**T* 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**T* 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**T*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]**T * 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]**T * 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**T * 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**T * 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**T*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