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