KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlatrs.f
Go to the documentation of this file.
1  SUBROUTINE dlatrs( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
2  $ CNORM, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * June 30, 1992
8 *
9 * .. Scalar Arguments ..
10  CHARACTER DIAG, NORMIN, TRANS, UPLO
11  INTEGER INFO, LDA, N
12  DOUBLE PRECISION SCALE
13 * ..
14 * .. Array Arguments ..
15  DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLATRS solves one of the triangular systems
22 *
23 * A *x = s*b or A'*x = s*b
24 *
25 * with scaling to prevent overflow. Here A is an upper or lower
26 * triangular matrix, A' denotes the transpose of A, x and b are
27 * n-element vectors, and s is a scaling factor, usually less than
28 * or equal to 1, chosen so that the components of x will be less than
29 * the overflow threshold. If the unscaled problem will not cause
30 * overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A
31 * is singular (A(j,j) = 0 for some j), then s is set to 0 and a
32 * non-trivial solution to A*x = 0 is returned.
33 *
34 * Arguments
35 * =========
36 *
37 * UPLO (input) CHARACTER*1
38 * Specifies whether the matrix A is upper or lower triangular.
39 * = 'U': Upper triangular
40 * = 'L': Lower triangular
41 *
42 * TRANS (input) CHARACTER*1
43 * Specifies the operation applied to A.
44 * = 'N': Solve A * x = s*b (No transpose)
45 * = 'T': Solve A'* x = s*b (Transpose)
46 * = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose)
47 *
48 * DIAG (input) CHARACTER*1
49 * Specifies whether or not the matrix A is unit triangular.
50 * = 'N': Non-unit triangular
51 * = 'U': Unit triangular
52 *
53 * NORMIN (input) CHARACTER*1
54 * Specifies whether CNORM has been set or not.
55 * = 'Y': CNORM contains the column norms on entry
56 * = 'N': CNORM is not set on entry. On exit, the norms will
57 * be computed and stored in CNORM.
58 *
59 * N (input) INTEGER
60 * The order of the matrix A. N >= 0.
61 *
62 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
63 * The triangular matrix A. If UPLO = 'U', the leading n by n
64 * upper triangular part of the array A contains the upper
65 * triangular matrix, and the strictly lower triangular part of
66 * A is not referenced. If UPLO = 'L', the leading n by n lower
67 * triangular part of the array A contains the lower triangular
68 * matrix, and the strictly upper triangular part of A is not
69 * referenced. If DIAG = 'U', the diagonal elements of A are
70 * also not referenced and are assumed to be 1.
71 *
72 * LDA (input) INTEGER
73 * The leading dimension of the array A. LDA >= max (1,N).
74 *
75 * X (input/output) DOUBLE PRECISION array, dimension (N)
76 * On entry, the right hand side b of the triangular system.
77 * On exit, X is overwritten by the solution vector x.
78 *
79 * SCALE (output) DOUBLE PRECISION
80 * The scaling factor s for the triangular system
81 * A * x = s*b or A'* x = s*b.
82 * If SCALE = 0, the matrix A is singular or badly scaled, and
83 * the vector x is an exact or approximate solution to A*x = 0.
84 *
85 * CNORM (input or output) DOUBLE PRECISION array, dimension (N)
86 *
87 * If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
88 * contains the norm of the off-diagonal part of the j-th column
89 * of A. If TRANS = 'N', CNORM(j) must be greater than or equal
90 * to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
91 * must be greater than or equal to the 1-norm.
92 *
93 * If NORMIN = 'N', CNORM is an output argument and CNORM(j)
94 * returns the 1-norm of the offdiagonal part of the j-th column
95 * of A.
96 *
97 * INFO (output) INTEGER
98 * = 0: successful exit
99 * < 0: if INFO = -k, the k-th argument had an illegal value
100 *
101 * Further Details
102 * ======= =======
103 *
104 * A rough bound on x is computed; if that is less than overflow, DTRSV
105 * is called, otherwise, specific code is used which checks for possible
106 * overflow or divide-by-zero at every operation.
107 *
108 * A columnwise scheme is used for solving A*x = b. The basic algorithm
109 * if A is lower triangular is
110 *
111 * x[1:n] := b[1:n]
112 * for j = 1, ..., n
113 * x(j) := x(j) / A(j,j)
114 * x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
115 * end
116 *
117 * Define bounds on the components of x after j iterations of the loop:
118 * M(j) = bound on x[1:j]
119 * G(j) = bound on x[j+1:n]
120 * Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
121 *
122 * Then for iteration j+1 we have
123 * M(j+1) <= G(j) / | A(j+1,j+1) |
124 * G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
125 * <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
126 *
127 * where CNORM(j+1) is greater than or equal to the infinity-norm of
128 * column j+1 of A, not counting the diagonal. Hence
129 *
130 * G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
131 * 1<=i<=j
132 * and
133 *
134 * |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
135 * 1<=i< j
136 *
137 * Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the
138 * reciprocal of the largest M(j), j=1,..,n, is larger than
139 * max(underflow, 1/overflow).
140 *
141 * The bound on x(j) is also used to determine when a step in the
142 * columnwise method can be performed without fear of overflow. If
143 * the computed bound is greater than a large constant, x is scaled to
144 * prevent overflow, but if the bound overflows, x is set to 0, x(j) to
145 * 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
146 *
147 * Similarly, a row-wise scheme is used to solve A'*x = b. The basic
148 * algorithm for A upper triangular is
149 *
150 * for j = 1, ..., n
151 * x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
152 * end
153 *
154 * We simultaneously compute two bounds
155 * G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
156 * M(j) = bound on x(i), 1<=i<=j
157 *
158 * The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
159 * add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
160 * Then the bound on x(j) is
161 *
162 * M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
163 *
164 * <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
165 * 1<=i<=j
166 *
167 * and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater
168 * than max(underflow, 1/overflow).
169 *
170 * =====================================================================
171 *
172 * .. Parameters ..
173  DOUBLE PRECISION ZERO, HALF, ONE
174  parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
175 * ..
176 * .. Local Scalars ..
177  LOGICAL NOTRAN, NOUNIT, UPPER
178  INTEGER I, IMAX, J, JFIRST, JINC, JLAST
179  DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
180  $ tmax, tscal, uscal, xbnd, xj, xmax
181 * ..
182 * .. External Functions ..
183  LOGICAL LSAME
184  INTEGER IDAMAX
185  DOUBLE PRECISION DASUM, DDOT, DLAMCH
186  EXTERNAL lsame, idamax, dasum, ddot, dlamch
187 * ..
188 * .. External Subroutines ..
189  EXTERNAL daxpy, dscal, dtrsv, xerbla
190 * ..
191 * .. Intrinsic Functions ..
192  INTRINSIC abs, max, min
193 * ..
194 * .. Executable Statements ..
195 *
196  info = 0
197  upper = lsame( uplo, 'U' )
198  notran = lsame( trans, 'N' )
199  nounit = lsame( diag, 'N' )
200 *
201 * Test the input parameters.
202 *
203  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
204  info = -1
205  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
206  $ lsame( trans, 'C' ) ) THEN
207  info = -2
208  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
209  info = -3
210  ELSE IF( .NOT.lsame( normin, 'Y' ) .AND. .NOT.
211  $ lsame( normin, 'N' ) ) THEN
212  info = -4
213  ELSE IF( n.LT.0 ) THEN
214  info = -5
215  ELSE IF( lda.LT.max( 1, n ) ) THEN
216  info = -7
217  END IF
218  IF( info.NE.0 ) THEN
219  CALL xerbla( 'DLATRS', -info )
220  RETURN
221  END IF
222 *
223 * Quick return if possible
224 *
225  IF( n.EQ.0 )
226  $ RETURN
227 *
228 * Determine machine dependent parameters to control overflow.
229 *
230  smlnum = dlamch( 'Safe minimum' ) / dlamch( 'Precision' )
231  bignum = one / smlnum
232  scale = one
233 *
234  IF( lsame( normin, 'N' ) ) THEN
235 *
236 * Compute the 1-norm of each column, not including the diagonal.
237 *
238  IF( upper ) THEN
239 *
240 * A is upper triangular.
241 *
242  DO 10 j = 1, n
243  cnorm( j ) = dasum( j-1, a( 1, j ), 1 )
244  10 CONTINUE
245  ELSE
246 *
247 * A is lower triangular.
248 *
249  DO 20 j = 1, n - 1
250  cnorm( j ) = dasum( n-j, a( j+1, j ), 1 )
251  20 CONTINUE
252  cnorm( n ) = zero
253  END IF
254  END IF
255 *
256 * Scale the column norms by TSCAL if the maximum element in CNORM is
257 * greater than BIGNUM.
258 *
259  imax = idamax( n, cnorm, 1 )
260  tmax = cnorm( imax )
261  IF( tmax.LE.bignum ) THEN
262  tscal = one
263  ELSE
264  tscal = one / ( smlnum*tmax )
265  CALL dscal( n, tscal, cnorm, 1 )
266  END IF
267 *
268 * Compute a bound on the computed solution vector to see if the
269 * Level 2 BLAS routine DTRSV can be used.
270 *
271  j = idamax( n, x, 1 )
272  xmax = abs( x( j ) )
273  xbnd = xmax
274  IF( notran ) THEN
275 *
276 * Compute the growth in A * x = b.
277 *
278  IF( upper ) THEN
279  jfirst = n
280  jlast = 1
281  jinc = -1
282  ELSE
283  jfirst = 1
284  jlast = n
285  jinc = 1
286  END IF
287 *
288  IF( tscal.NE.one ) THEN
289  grow = zero
290  GO TO 50
291  END IF
292 *
293  IF( nounit ) THEN
294 *
295 * A is non-unit triangular.
296 *
297 * Compute GROW = 1/G(j) and XBND = 1/M(j).
298 * Initially, G(0) = max{x(i), i=1,...,n}.
299 *
300  grow = one / max( xbnd, smlnum )
301  xbnd = grow
302  DO 30 j = jfirst, jlast, jinc
303 *
304 * Exit the loop if the growth factor is too small.
305 *
306  IF( grow.LE.smlnum )
307  $ GO TO 50
308 *
309 * M(j) = G(j-1) / abs(A(j,j))
310 *
311  tjj = abs( a( j, j ) )
312  xbnd = min( xbnd, min( one, tjj )*grow )
313  IF( tjj+cnorm( j ).GE.smlnum ) THEN
314 *
315 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
316 *
317  grow = grow*( tjj / ( tjj+cnorm( j ) ) )
318  ELSE
319 *
320 * G(j) could overflow, set GROW to 0.
321 *
322  grow = zero
323  END IF
324  30 CONTINUE
325  grow = xbnd
326  ELSE
327 *
328 * A is unit triangular.
329 *
330 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
331 *
332  grow = min( one, one / max( xbnd, smlnum ) )
333  DO 40 j = jfirst, jlast, jinc
334 *
335 * Exit the loop if the growth factor is too small.
336 *
337  IF( grow.LE.smlnum )
338  $ GO TO 50
339 *
340 * G(j) = G(j-1)*( 1 + CNORM(j) )
341 *
342  grow = grow*( one / ( one+cnorm( j ) ) )
343  40 CONTINUE
344  END IF
345  50 CONTINUE
346 *
347  ELSE
348 *
349 * Compute the growth in A' * x = b.
350 *
351  IF( upper ) THEN
352  jfirst = 1
353  jlast = n
354  jinc = 1
355  ELSE
356  jfirst = n
357  jlast = 1
358  jinc = -1
359  END IF
360 *
361  IF( tscal.NE.one ) THEN
362  grow = zero
363  GO TO 80
364  END IF
365 *
366  IF( nounit ) THEN
367 *
368 * A is non-unit triangular.
369 *
370 * Compute GROW = 1/G(j) and XBND = 1/M(j).
371 * Initially, M(0) = max{x(i), i=1,...,n}.
372 *
373  grow = one / max( xbnd, smlnum )
374  xbnd = grow
375  DO 60 j = jfirst, jlast, jinc
376 *
377 * Exit the loop if the growth factor is too small.
378 *
379  IF( grow.LE.smlnum )
380  $ GO TO 80
381 *
382 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
383 *
384  xj = one + cnorm( j )
385  grow = min( grow, xbnd / xj )
386 *
387 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
388 *
389  tjj = abs( a( j, j ) )
390  IF( xj.GT.tjj )
391  $ xbnd = xbnd*( tjj / xj )
392  60 CONTINUE
393  grow = min( grow, xbnd )
394  ELSE
395 *
396 * A is unit triangular.
397 *
398 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
399 *
400  grow = min( one, one / max( xbnd, smlnum ) )
401  DO 70 j = jfirst, jlast, jinc
402 *
403 * Exit the loop if the growth factor is too small.
404 *
405  IF( grow.LE.smlnum )
406  $ GO TO 80
407 *
408 * G(j) = ( 1 + CNORM(j) )*G(j-1)
409 *
410  xj = one + cnorm( j )
411  grow = grow / xj
412  70 CONTINUE
413  END IF
414  80 CONTINUE
415  END IF
416 *
417  IF( ( grow*tscal ).GT.smlnum ) THEN
418 *
419 * Use the Level 2 BLAS solve if the reciprocal of the bound on
420 * elements of X is not too small.
421 *
422  CALL dtrsv( uplo, trans, diag, n, a, lda, x, 1 )
423  ELSE
424 *
425 * Use a Level 1 BLAS solve, scaling intermediate results.
426 *
427  IF( xmax.GT.bignum ) THEN
428 *
429 * Scale X so that its components are less than or equal to
430 * BIGNUM in absolute value.
431 *
432  scale = bignum / xmax
433  CALL dscal( n, scale, x, 1 )
434  xmax = bignum
435  END IF
436 *
437  IF( notran ) THEN
438 *
439 * Solve A * x = b
440 *
441  DO 110 j = jfirst, jlast, jinc
442 *
443 * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
444 *
445  xj = abs( x( j ) )
446  IF( nounit ) THEN
447  tjjs = a( j, j )*tscal
448  ELSE
449  tjjs = tscal
450  IF( tscal.EQ.one )
451  $ GO TO 100
452  END IF
453  tjj = abs( tjjs )
454  IF( tjj.GT.smlnum ) THEN
455 *
456 * abs(A(j,j)) > SMLNUM:
457 *
458  IF( tjj.LT.one ) THEN
459  IF( xj.GT.tjj*bignum ) THEN
460 *
461 * Scale x by 1/b(j).
462 *
463  rec = one / xj
464  CALL dscal( n, rec, x, 1 )
465  scale = scale*rec
466  xmax = xmax*rec
467  END IF
468  END IF
469  x( j ) = x( j ) / tjjs
470  xj = abs( x( j ) )
471  ELSE IF( tjj.GT.zero ) THEN
472 *
473 * 0 < abs(A(j,j)) <= SMLNUM:
474 *
475  IF( xj.GT.tjj*bignum ) THEN
476 *
477 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
478 * to avoid overflow when dividing by A(j,j).
479 *
480  rec = ( tjj*bignum ) / xj
481  IF( cnorm( j ).GT.one ) THEN
482 *
483 * Scale by 1/CNORM(j) to avoid overflow when
484 * multiplying x(j) times column j.
485 *
486  rec = rec / cnorm( j )
487  END IF
488  CALL dscal( n, rec, x, 1 )
489  scale = scale*rec
490  xmax = xmax*rec
491  END IF
492  x( j ) = x( j ) / tjjs
493  xj = abs( x( j ) )
494  ELSE
495 *
496 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
497 * scale = 0, and compute a solution to A*x = 0.
498 *
499  DO 90 i = 1, n
500  x( i ) = zero
501  90 CONTINUE
502  x( j ) = one
503  xj = one
504  scale = zero
505  xmax = zero
506  END IF
507  100 CONTINUE
508 *
509 * Scale x if necessary to avoid overflow when adding a
510 * multiple of column j of A.
511 *
512  IF( xj.GT.one ) THEN
513  rec = one / xj
514  IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
515 *
516 * Scale x by 1/(2*abs(x(j))).
517 *
518  rec = rec*half
519  CALL dscal( n, rec, x, 1 )
520  scale = scale*rec
521  END IF
522  ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
523 *
524 * Scale x by 1/2.
525 *
526  CALL dscal( n, half, x, 1 )
527  scale = scale*half
528  END IF
529 *
530  IF( upper ) THEN
531  IF( j.GT.1 ) THEN
532 *
533 * Compute the update
534 * x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
535 *
536  CALL daxpy( j-1, -x( j )*tscal, a( 1, j ), 1, x,
537  $ 1 )
538  i = idamax( j-1, x, 1 )
539  xmax = abs( x( i ) )
540  END IF
541  ELSE
542  IF( j.LT.n ) THEN
543 *
544 * Compute the update
545 * x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
546 *
547  CALL daxpy( n-j, -x( j )*tscal, a( j+1, j ), 1,
548  $ x( j+1 ), 1 )
549  i = j + idamax( n-j, x( j+1 ), 1 )
550  xmax = abs( x( i ) )
551  END IF
552  END IF
553  110 CONTINUE
554 *
555  ELSE
556 *
557 * Solve A' * x = b
558 *
559  DO 160 j = jfirst, jlast, jinc
560 *
561 * Compute x(j) = b(j) - sum A(k,j)*x(k).
562 * k<>j
563 *
564  xj = abs( x( j ) )
565  uscal = tscal
566  rec = one / max( xmax, one )
567  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
568 *
569 * If x(j) could overflow, scale x by 1/(2*XMAX).
570 *
571  rec = rec*half
572  IF( nounit ) THEN
573  tjjs = a( j, j )*tscal
574  ELSE
575  tjjs = tscal
576  END IF
577  tjj = abs( tjjs )
578  IF( tjj.GT.one ) THEN
579 *
580 * Divide by A(j,j) when scaling x if A(j,j) > 1.
581 *
582  rec = min( one, rec*tjj )
583  uscal = uscal / tjjs
584  END IF
585  IF( rec.LT.one ) THEN
586  CALL dscal( n, rec, x, 1 )
587  scale = scale*rec
588  xmax = xmax*rec
589  END IF
590  END IF
591 *
592  sumj = zero
593  IF( uscal.EQ.one ) THEN
594 *
595 * If the scaling needed for A in the dot product is 1,
596 * call DDOT to perform the dot product.
597 *
598  IF( upper ) THEN
599  sumj = ddot( j-1, a( 1, j ), 1, x, 1 )
600  ELSE IF( j.LT.n ) THEN
601  sumj = ddot( n-j, a( j+1, j ), 1, x( j+1 ), 1 )
602  END IF
603  ELSE
604 *
605 * Otherwise, use in-line code for the dot product.
606 *
607  IF( upper ) THEN
608  DO 120 i = 1, j - 1
609  sumj = sumj + ( a( i, j )*uscal )*x( i )
610  120 CONTINUE
611  ELSE IF( j.LT.n ) THEN
612  DO 130 i = j + 1, n
613  sumj = sumj + ( a( i, j )*uscal )*x( i )
614  130 CONTINUE
615  END IF
616  END IF
617 *
618  IF( uscal.EQ.tscal ) THEN
619 *
620 * Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
621 * was not used to scale the dotproduct.
622 *
623  x( j ) = x( j ) - sumj
624  xj = abs( x( j ) )
625  IF( nounit ) THEN
626  tjjs = a( j, j )*tscal
627  ELSE
628  tjjs = tscal
629  IF( tscal.EQ.one )
630  $ GO TO 150
631  END IF
632 *
633 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
634 *
635  tjj = abs( tjjs )
636  IF( tjj.GT.smlnum ) THEN
637 *
638 * abs(A(j,j)) > SMLNUM:
639 *
640  IF( tjj.LT.one ) THEN
641  IF( xj.GT.tjj*bignum ) THEN
642 *
643 * Scale X by 1/abs(x(j)).
644 *
645  rec = one / xj
646  CALL dscal( n, rec, x, 1 )
647  scale = scale*rec
648  xmax = xmax*rec
649  END IF
650  END IF
651  x( j ) = x( j ) / tjjs
652  ELSE IF( tjj.GT.zero ) THEN
653 *
654 * 0 < abs(A(j,j)) <= SMLNUM:
655 *
656  IF( xj.GT.tjj*bignum ) THEN
657 *
658 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
659 *
660  rec = ( tjj*bignum ) / xj
661  CALL dscal( n, rec, x, 1 )
662  scale = scale*rec
663  xmax = xmax*rec
664  END IF
665  x( j ) = x( j ) / tjjs
666  ELSE
667 *
668 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
669 * scale = 0, and compute a solution to A'*x = 0.
670 *
671  DO 140 i = 1, n
672  x( i ) = zero
673  140 CONTINUE
674  x( j ) = one
675  scale = zero
676  xmax = zero
677  END IF
678  150 CONTINUE
679  ELSE
680 *
681 * Compute x(j) := x(j) / A(j,j) - sumj if the dot
682 * product has already been divided by 1/A(j,j).
683 *
684  x( j ) = x( j ) / tjjs - sumj
685  END IF
686  xmax = max( xmax, abs( x( j ) ) )
687  160 CONTINUE
688  END IF
689  scale = scale / tscal
690  END IF
691 *
692 * Scale the column norms by 1/TSCAL for return.
693 *
694  IF( tscal.NE.one ) THEN
695  CALL dscal( n, one / tscal, cnorm, 1 )
696  END IF
697 *
698  RETURN
699 *
700 * End of DLATRS
701 *
702  END
subroutine daxpy(n, da, dx, incx, dy, incy)
Definition: daxpy.f:2
subroutine dlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
Definition: dlatrs.f:3
subroutine dscal(n, da, dx, incx)
Definition: dscal.f:2
subroutine dtrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
Definition: dtrsv.f:2
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2