LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dlatps.f
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
dlatps
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
dtpsv
subroutine dtpsv(UPLO, TRANS, DIAG, N, AP, X, INCX)
DTPSV
Definition: dtpsv.f:144
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79