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