LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlatps.f
Go to the documentation of this file.
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 latps
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
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dlatps(uplo, trans, diag, normin, n, ap, x, scale, cnorm, info)
DLATPS solves a triangular system of equations with the matrix held in packed storage.
Definition dlatps.f:229
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dtpsv(uplo, trans, diag, n, ap, x, incx)
DTPSV
Definition dtpsv.f:144