LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
clatrs.f
1 *> \brief \b CLATRS 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 CLATRS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatrs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatrs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatrs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLATRS( 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 * REAL SCALE
28 * ..
29 * .. Array Arguments ..
30 * REAL CNORM( * )
31 * COMPLEX A( LDA, * ), X( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CLATRS 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 *> CTRSV 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 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 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 REAL
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 REAL 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 latrs
163 *
164 *> \par Further Details:
165 * =====================
166 *>
167 *> \verbatim
168 *>
169 *> A rough bound on x is computed; if that is less than overflow, CTRSV
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 CTRSV 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 CTRSV if 1/M(n) and 1/G(n) are both greater
233 *> than max(underflow, 1/overflow).
234 *> \endverbatim
235 *>
236 * =====================================================================
237  SUBROUTINE clatrs( 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  REAL SCALE
248 * ..
249 * .. Array Arguments ..
250  REAL CNORM( * )
251  COMPLEX A( lda, * ), X( * )
252 * ..
253 *
254 * =====================================================================
255 *
256 * .. Parameters ..
257  REAL ZERO, HALF, ONE, TWO
258  parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
259  $ two = 2.0e+0 )
260 * ..
261 * .. Local Scalars ..
262  LOGICAL NOTRAN, NOUNIT, UPPER
263  INTEGER I, IMAX, J, JFIRST, JINC, JLAST
264  REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
265  $ xbnd, xj, xmax
266  COMPLEX CSUMJ, TJJS, USCAL, ZDUM
267 * ..
268 * .. External Functions ..
269  LOGICAL LSAME
270  INTEGER ICAMAX, ISAMAX
271  REAL SCASUM, SLAMCH
272  COMPLEX CDOTC, CDOTU, CLADIV
273  EXTERNAL lsame, icamax, isamax, scasum, slamch, cdotc,
274  $ cdotu, cladiv
275 * ..
276 * .. External Subroutines ..
277  EXTERNAL caxpy, csscal, ctrsv, sscal, xerbla
278 * ..
279 * .. Intrinsic Functions ..
280  INTRINSIC abs, aimag, cmplx, conjg, max, min, real
281 * ..
282 * .. Statement Functions ..
283  REAL CABS1, CABS2
284 * ..
285 * .. Statement Function definitions ..
286  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( AIMAG( zdum ) )
287  cabs2( zdum ) = abs( REAL( ZDUM ) / 2. ) +
288  $ abs( aimag( zdum ) / 2. )
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( 'CLATRS', -info )
316  RETURN
317  END IF
318 *
319 * Quick return if possible
320 *
321  scale = one
322  IF( n.EQ.0 )
323  $ RETURN
324 *
325 * Determine machine dependent parameters to control overflow.
326 *
327  smlnum = slamch( 'Safe minimum' ) / slamch( 'Precision' )
328  bignum = one / smlnum
329 *
330  IF( lsame( normin, 'N' ) ) THEN
331 *
332 * Compute the 1-norm of each column, not including the diagonal.
333 *
334  IF( upper ) THEN
335 *
336 * A is upper triangular.
337 *
338  DO 10 j = 1, n
339  cnorm( j ) = scasum( j-1, a( 1, j ), 1 )
340  10 CONTINUE
341  ELSE
342 *
343 * A is lower triangular.
344 *
345  DO 20 j = 1, n - 1
346  cnorm( j ) = scasum( n-j, a( j+1, j ), 1 )
347  20 CONTINUE
348  cnorm( n ) = zero
349  END IF
350  END IF
351 *
352 * Scale the column norms by TSCAL if the maximum element in CNORM is
353 * greater than BIGNUM/2.
354 *
355  imax = isamax( n, cnorm, 1 )
356  tmax = cnorm( imax )
357  IF( tmax.LE.bignum*half ) THEN
358  tscal = one
359  ELSE
360 *
361 * Avoid NaN generation if entries in CNORM exceed the
362 * overflow threshold
363 *
364  IF ( tmax.LE.slamch('Overflow') ) THEN
365 * Case 1: All entries in CNORM are valid floating-point numbers
366  tscal = half / ( smlnum*tmax )
367  CALL sscal( n, tscal, cnorm, 1 )
368  ELSE
369 * Case 2: At least one column norm of A cannot be
370 * represented as a floating-point number. Find the
371 * maximum offdiagonal absolute value
372 * max( |Re(A(I,J))|, |Im(A(I,J)| ). If this entry is
373 * not +/- Infinity, use this value as TSCAL.
374  tmax = zero
375  IF( upper ) THEN
376 *
377 * A is upper triangular.
378 *
379  DO j = 2, n
380  DO i = 1, j - 1
381  tmax = max( tmax, abs( REAL( A( I, J ) ) ),
382  $ abs( aimag(a( i, j ) ) ) )
383  END DO
384  END DO
385  ELSE
386 *
387 * A is lower triangular.
388 *
389  DO j = 1, n - 1
390  DO i = j + 1, n
391  tmax = max( tmax, abs( REAL( A( I, J ) ) ),
392  $ abs( aimag(a( i, j ) ) ) )
393  END DO
394  END DO
395  END IF
396 *
397  IF( tmax.LE.slamch('Overflow') ) THEN
398  tscal = one / ( smlnum*tmax )
399  DO j = 1, n
400  IF( cnorm( j ).LE.slamch('Overflow') ) THEN
401  cnorm( j ) = cnorm( j )*tscal
402  ELSE
403 * Recompute the 1-norm of each column without
404 * introducing Infinity in the summation.
405  tscal = two * tscal
406  cnorm( j ) = zero
407  IF( upper ) THEN
408  DO i = 1, j - 1
409  cnorm( j ) = cnorm( j ) +
410  $ tscal * cabs2( a( i, j ) )
411  END DO
412  ELSE
413  DO i = j + 1, n
414  cnorm( j ) = cnorm( j ) +
415  $ tscal * cabs2( a( i, j ) )
416  END DO
417  END IF
418  tscal = tscal * half
419  END IF
420  END DO
421  ELSE
422 * At least one entry of A is not a valid floating-point
423 * entry. Rely on TRSV to propagate Inf and NaN.
424  CALL ctrsv( uplo, trans, diag, n, a, lda, x, 1 )
425  RETURN
426  END IF
427  END IF
428  END IF
429 *
430 * Compute a bound on the computed solution vector to see if the
431 * Level 2 BLAS routine CTRSV can be used.
432 *
433  xmax = zero
434  DO 30 j = 1, n
435  xmax = max( xmax, cabs2( x( j ) ) )
436  30 CONTINUE
437  xbnd = xmax
438 *
439  IF( notran ) THEN
440 *
441 * Compute the growth in A * x = b.
442 *
443  IF( upper ) THEN
444  jfirst = n
445  jlast = 1
446  jinc = -1
447  ELSE
448  jfirst = 1
449  jlast = n
450  jinc = 1
451  END IF
452 *
453  IF( tscal.NE.one ) THEN
454  grow = zero
455  GO TO 60
456  END IF
457 *
458  IF( nounit ) THEN
459 *
460 * A is non-unit triangular.
461 *
462 * Compute GROW = 1/G(j) and XBND = 1/M(j).
463 * Initially, G(0) = max{x(i), i=1,...,n}.
464 *
465  grow = half / max( xbnd, smlnum )
466  xbnd = grow
467  DO 40 j = jfirst, jlast, jinc
468 *
469 * Exit the loop if the growth factor is too small.
470 *
471  IF( grow.LE.smlnum )
472  $ GO TO 60
473 *
474  tjjs = a( j, j )
475  tjj = cabs1( tjjs )
476 *
477  IF( tjj.GE.smlnum ) THEN
478 *
479 * M(j) = G(j-1) / abs(A(j,j))
480 *
481  xbnd = min( xbnd, min( one, tjj )*grow )
482  ELSE
483 *
484 * M(j) could overflow, set XBND to 0.
485 *
486  xbnd = zero
487  END IF
488 *
489  IF( tjj+cnorm( j ).GE.smlnum ) THEN
490 *
491 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
492 *
493  grow = grow*( tjj / ( tjj+cnorm( j ) ) )
494  ELSE
495 *
496 * G(j) could overflow, set GROW to 0.
497 *
498  grow = zero
499  END IF
500  40 CONTINUE
501  grow = xbnd
502  ELSE
503 *
504 * A is unit triangular.
505 *
506 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
507 *
508  grow = min( one, half / max( xbnd, smlnum ) )
509  DO 50 j = jfirst, jlast, jinc
510 *
511 * Exit the loop if the growth factor is too small.
512 *
513  IF( grow.LE.smlnum )
514  $ GO TO 60
515 *
516 * G(j) = G(j-1)*( 1 + CNORM(j) )
517 *
518  grow = grow*( one / ( one+cnorm( j ) ) )
519  50 CONTINUE
520  END IF
521  60 CONTINUE
522 *
523  ELSE
524 *
525 * Compute the growth in A**T * x = b or A**H * x = b.
526 *
527  IF( upper ) THEN
528  jfirst = 1
529  jlast = n
530  jinc = 1
531  ELSE
532  jfirst = n
533  jlast = 1
534  jinc = -1
535  END IF
536 *
537  IF( tscal.NE.one ) THEN
538  grow = zero
539  GO TO 90
540  END IF
541 *
542  IF( nounit ) THEN
543 *
544 * A is non-unit triangular.
545 *
546 * Compute GROW = 1/G(j) and XBND = 1/M(j).
547 * Initially, M(0) = max{x(i), i=1,...,n}.
548 *
549  grow = half / max( xbnd, smlnum )
550  xbnd = grow
551  DO 70 j = jfirst, jlast, jinc
552 *
553 * Exit the loop if the growth factor is too small.
554 *
555  IF( grow.LE.smlnum )
556  $ GO TO 90
557 *
558 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
559 *
560  xj = one + cnorm( j )
561  grow = min( grow, xbnd / xj )
562 *
563  tjjs = a( j, j )
564  tjj = cabs1( tjjs )
565 *
566  IF( tjj.GE.smlnum ) THEN
567 *
568 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
569 *
570  IF( xj.GT.tjj )
571  $ xbnd = xbnd*( tjj / xj )
572  ELSE
573 *
574 * M(j) could overflow, set XBND to 0.
575 *
576  xbnd = zero
577  END IF
578  70 CONTINUE
579  grow = min( grow, xbnd )
580  ELSE
581 *
582 * A is unit triangular.
583 *
584 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
585 *
586  grow = min( one, half / max( xbnd, smlnum ) )
587  DO 80 j = jfirst, jlast, jinc
588 *
589 * Exit the loop if the growth factor is too small.
590 *
591  IF( grow.LE.smlnum )
592  $ GO TO 90
593 *
594 * G(j) = ( 1 + CNORM(j) )*G(j-1)
595 *
596  xj = one + cnorm( j )
597  grow = grow / xj
598  80 CONTINUE
599  END IF
600  90 CONTINUE
601  END IF
602 *
603  IF( ( grow*tscal ).GT.smlnum ) THEN
604 *
605 * Use the Level 2 BLAS solve if the reciprocal of the bound on
606 * elements of X is not too small.
607 *
608  CALL ctrsv( uplo, trans, diag, n, a, lda, x, 1 )
609  ELSE
610 *
611 * Use a Level 1 BLAS solve, scaling intermediate results.
612 *
613  IF( xmax.GT.bignum*half ) THEN
614 *
615 * Scale X so that its components are less than or equal to
616 * BIGNUM in absolute value.
617 *
618  scale = ( bignum*half ) / xmax
619  CALL csscal( n, scale, x, 1 )
620  xmax = bignum
621  ELSE
622  xmax = xmax*two
623  END IF
624 *
625  IF( notran ) THEN
626 *
627 * Solve A * x = b
628 *
629  DO 110 j = jfirst, jlast, jinc
630 *
631 * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
632 *
633  xj = cabs1( x( j ) )
634  IF( nounit ) THEN
635  tjjs = a( j, j )*tscal
636  ELSE
637  tjjs = tscal
638  IF( tscal.EQ.one )
639  $ GO TO 105
640  END IF
641  tjj = cabs1( tjjs )
642  IF( tjj.GT.smlnum ) THEN
643 *
644 * abs(A(j,j)) > SMLNUM:
645 *
646  IF( tjj.LT.one ) THEN
647  IF( xj.GT.tjj*bignum ) THEN
648 *
649 * Scale x by 1/b(j).
650 *
651  rec = one / xj
652  CALL csscal( n, rec, x, 1 )
653  scale = scale*rec
654  xmax = xmax*rec
655  END IF
656  END IF
657  x( j ) = cladiv( x( j ), tjjs )
658  xj = cabs1( x( j ) )
659  ELSE IF( tjj.GT.zero ) THEN
660 *
661 * 0 < abs(A(j,j)) <= SMLNUM:
662 *
663  IF( xj.GT.tjj*bignum ) THEN
664 *
665 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
666 * to avoid overflow when dividing by A(j,j).
667 *
668  rec = ( tjj*bignum ) / xj
669  IF( cnorm( j ).GT.one ) THEN
670 *
671 * Scale by 1/CNORM(j) to avoid overflow when
672 * multiplying x(j) times column j.
673 *
674  rec = rec / cnorm( j )
675  END IF
676  CALL csscal( n, rec, x, 1 )
677  scale = scale*rec
678  xmax = xmax*rec
679  END IF
680  x( j ) = cladiv( x( j ), tjjs )
681  xj = cabs1( x( j ) )
682  ELSE
683 *
684 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
685 * scale = 0, and compute a solution to A*x = 0.
686 *
687  DO 100 i = 1, n
688  x( i ) = zero
689  100 CONTINUE
690  x( j ) = one
691  xj = one
692  scale = zero
693  xmax = zero
694  END IF
695  105 CONTINUE
696 *
697 * Scale x if necessary to avoid overflow when adding a
698 * multiple of column j of A.
699 *
700  IF( xj.GT.one ) THEN
701  rec = one / xj
702  IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
703 *
704 * Scale x by 1/(2*abs(x(j))).
705 *
706  rec = rec*half
707  CALL csscal( n, rec, x, 1 )
708  scale = scale*rec
709  END IF
710  ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
711 *
712 * Scale x by 1/2.
713 *
714  CALL csscal( n, half, x, 1 )
715  scale = scale*half
716  END IF
717 *
718  IF( upper ) THEN
719  IF( j.GT.1 ) THEN
720 *
721 * Compute the update
722 * x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
723 *
724  CALL caxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
725  $ 1 )
726  i = icamax( j-1, x, 1 )
727  xmax = cabs1( x( i ) )
728  END IF
729  ELSE
730  IF( j.LT.n ) THEN
731 *
732 * Compute the update
733 * x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
734 *
735  CALL caxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
736  $ x( j+1 ), 1 )
737  i = j + icamax( n-j, x( j+1 ), 1 )
738  xmax = cabs1( x( i ) )
739  END IF
740  END IF
741  110 CONTINUE
742 *
743  ELSE IF( lsame( trans, 'T' ) ) THEN
744 *
745 * Solve A**T * x = b
746 *
747  DO 150 j = jfirst, jlast, jinc
748 *
749 * Compute x(j) = b(j) - sum A(k,j)*x(k).
750 * k<>j
751 *
752  xj = cabs1( x( j ) )
753  uscal = tscal
754  rec = one / max( xmax, one )
755  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
756 *
757 * If x(j) could overflow, scale x by 1/(2*XMAX).
758 *
759  rec = rec*half
760  IF( nounit ) THEN
761  tjjs = a( j, j )*tscal
762  ELSE
763  tjjs = tscal
764  END IF
765  tjj = cabs1( tjjs )
766  IF( tjj.GT.one ) THEN
767 *
768 * Divide by A(j,j) when scaling x if A(j,j) > 1.
769 *
770  rec = min( one, rec*tjj )
771  uscal = cladiv( uscal, tjjs )
772  END IF
773  IF( rec.LT.one ) THEN
774  CALL csscal( n, rec, x, 1 )
775  scale = scale*rec
776  xmax = xmax*rec
777  END IF
778  END IF
779 *
780  csumj = zero
781  IF( uscal.EQ.cmplx( one ) ) THEN
782 *
783 * If the scaling needed for A in the dot product is 1,
784 * call CDOTU to perform the dot product.
785 *
786  IF( upper ) THEN
787  csumj = cdotu( j-1, a( 1, j ), 1, x, 1 )
788  ELSE IF( j.LT.n ) THEN
789  csumj = cdotu( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
790  END IF
791  ELSE
792 *
793 * Otherwise, use in-line code for the dot product.
794 *
795  IF( upper ) THEN
796  DO 120 i = 1, j - 1
797  csumj = csumj + ( a( i, j )*uscal )*x( i )
798  120 CONTINUE
799  ELSE IF( j.LT.n ) THEN
800  DO 130 i = j + 1, n
801  csumj = csumj + ( a( i, j )*uscal )*x( i )
802  130 CONTINUE
803  END IF
804  END IF
805 *
806  IF( uscal.EQ.cmplx( tscal ) ) THEN
807 *
808 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
809 * was not used to scale the dotproduct.
810 *
811  x( j ) = x( j ) - csumj
812  xj = cabs1( x( j ) )
813  IF( nounit ) THEN
814  tjjs = a( j, j )*tscal
815  ELSE
816  tjjs = tscal
817  IF( tscal.EQ.one )
818  $ GO TO 145
819  END IF
820 *
821 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
822 *
823  tjj = cabs1( tjjs )
824  IF( tjj.GT.smlnum ) THEN
825 *
826 * abs(A(j,j)) > SMLNUM:
827 *
828  IF( tjj.LT.one ) THEN
829  IF( xj.GT.tjj*bignum ) THEN
830 *
831 * Scale X by 1/abs(x(j)).
832 *
833  rec = one / xj
834  CALL csscal( n, rec, x, 1 )
835  scale = scale*rec
836  xmax = xmax*rec
837  END IF
838  END IF
839  x( j ) = cladiv( x( j ), tjjs )
840  ELSE IF( tjj.GT.zero ) THEN
841 *
842 * 0 < abs(A(j,j)) <= SMLNUM:
843 *
844  IF( xj.GT.tjj*bignum ) THEN
845 *
846 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
847 *
848  rec = ( tjj*bignum ) / xj
849  CALL csscal( n, rec, x, 1 )
850  scale = scale*rec
851  xmax = xmax*rec
852  END IF
853  x( j ) = cladiv( x( j ), tjjs )
854  ELSE
855 *
856 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
857 * scale = 0 and compute a solution to A**T *x = 0.
858 *
859  DO 140 i = 1, n
860  x( i ) = zero
861  140 CONTINUE
862  x( j ) = one
863  scale = zero
864  xmax = zero
865  END IF
866  145 CONTINUE
867  ELSE
868 *
869 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
870 * product has already been divided by 1/A(j,j).
871 *
872  x( j ) = cladiv( x( j ), tjjs ) - csumj
873  END IF
874  xmax = max( xmax, cabs1( x( j ) ) )
875  150 CONTINUE
876 *
877  ELSE
878 *
879 * Solve A**H * x = b
880 *
881  DO 190 j = jfirst, jlast, jinc
882 *
883 * Compute x(j) = b(j) - sum A(k,j)*x(k).
884 * k<>j
885 *
886  xj = cabs1( x( j ) )
887  uscal = tscal
888  rec = one / max( xmax, one )
889  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
890 *
891 * If x(j) could overflow, scale x by 1/(2*XMAX).
892 *
893  rec = rec*half
894  IF( nounit ) THEN
895  tjjs = conjg( a( j, j ) )*tscal
896  ELSE
897  tjjs = tscal
898  END IF
899  tjj = cabs1( tjjs )
900  IF( tjj.GT.one ) THEN
901 *
902 * Divide by A(j,j) when scaling x if A(j,j) > 1.
903 *
904  rec = min( one, rec*tjj )
905  uscal = cladiv( uscal, tjjs )
906  END IF
907  IF( rec.LT.one ) THEN
908  CALL csscal( n, rec, x, 1 )
909  scale = scale*rec
910  xmax = xmax*rec
911  END IF
912  END IF
913 *
914  csumj = zero
915  IF( uscal.EQ.cmplx( one ) ) THEN
916 *
917 * If the scaling needed for A in the dot product is 1,
918 * call CDOTC to perform the dot product.
919 *
920  IF( upper ) THEN
921  csumj = cdotc( j-1, a( 1, j ), 1, x, 1 )
922  ELSE IF( j.LT.n ) THEN
923  csumj = cdotc( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
924  END IF
925  ELSE
926 *
927 * Otherwise, use in-line code for the dot product.
928 *
929  IF( upper ) THEN
930  DO 160 i = 1, j - 1
931  csumj = csumj + ( conjg( a( i, j ) )*uscal )*
932  $ x( i )
933  160 CONTINUE
934  ELSE IF( j.LT.n ) THEN
935  DO 170 i = j + 1, n
936  csumj = csumj + ( conjg( a( i, j ) )*uscal )*
937  $ x( i )
938  170 CONTINUE
939  END IF
940  END IF
941 *
942  IF( uscal.EQ.cmplx( tscal ) ) THEN
943 *
944 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
945 * was not used to scale the dotproduct.
946 *
947  x( j ) = x( j ) - csumj
948  xj = cabs1( x( j ) )
949  IF( nounit ) THEN
950  tjjs = conjg( a( j, j ) )*tscal
951  ELSE
952  tjjs = tscal
953  IF( tscal.EQ.one )
954  $ GO TO 185
955  END IF
956 *
957 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
958 *
959  tjj = cabs1( tjjs )
960  IF( tjj.GT.smlnum ) THEN
961 *
962 * abs(A(j,j)) > SMLNUM:
963 *
964  IF( tjj.LT.one ) THEN
965  IF( xj.GT.tjj*bignum ) THEN
966 *
967 * Scale X by 1/abs(x(j)).
968 *
969  rec = one / xj
970  CALL csscal( n, rec, x, 1 )
971  scale = scale*rec
972  xmax = xmax*rec
973  END IF
974  END IF
975  x( j ) = cladiv( x( j ), tjjs )
976  ELSE IF( tjj.GT.zero ) THEN
977 *
978 * 0 < abs(A(j,j)) <= SMLNUM:
979 *
980  IF( xj.GT.tjj*bignum ) THEN
981 *
982 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
983 *
984  rec = ( tjj*bignum ) / xj
985  CALL csscal( n, rec, x, 1 )
986  scale = scale*rec
987  xmax = xmax*rec
988  END IF
989  x( j ) = cladiv( x( j ), tjjs )
990  ELSE
991 *
992 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
993 * scale = 0 and compute a solution to A**H *x = 0.
994 *
995  DO 180 i = 1, n
996  x( i ) = zero
997  180 CONTINUE
998  x( j ) = one
999  scale = zero
1000  xmax = zero
1001  END IF
1002  185 CONTINUE
1003  ELSE
1004 *
1005 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
1006 * product has already been divided by 1/A(j,j).
1007 *
1008  x( j ) = cladiv( x( j ), tjjs ) - csumj
1009  END IF
1010  xmax = max( xmax, cabs1( x( j ) ) )
1011  190 CONTINUE
1012  END IF
1013  scale = scale / tscal
1014  END IF
1015 *
1016 * Scale the column norms by 1/TSCAL for return.
1017 *
1018  IF( tscal.NE.one ) THEN
1019  CALL sscal( n, one / tscal, cnorm, 1 )
1020  END IF
1021 *
1022  RETURN
1023 *
1024 * End of CLATRS
1025 *
1026  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
ctrsv
subroutine ctrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRSV
Definition: ctrsv.f:149
clatrs
subroutine clatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
CLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
Definition: clatrs.f:239
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88