ergo
template_lapack_latrs.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /* This file belongs to the template_lapack part of the Ergo source
31 * code. The source files in the template_lapack directory are modified
32 * versions of files originally distributed as CLAPACK, see the
33 * Copyright/license notice in the file template_lapack/COPYING.
34 */
35
36
37#ifndef TEMPLATE_LAPACK_LATRS_HEADER
38#define TEMPLATE_LAPACK_LATRS_HEADER
39
40
41template<class Treal>
42int template_lapack_latrs(const char *uplo, const char *trans, const char *diag, const char *
43 normin, const integer *n, const Treal *a, const integer *lda, Treal *x,
44 Treal *scale, Treal *cnorm, integer *info)
45{
46/* -- LAPACK auxiliary routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 June 30, 1992
50
51
52 Purpose
53 =======
54
55 DLATRS solves one of the triangular systems
56
57 A *x = s*b or A'*x = s*b
58
59 with scaling to prevent overflow. Here A is an upper or lower
60 triangular matrix, A' denotes the transpose of A, x and b are
61 n-element vectors, and s is a scaling factor, usually less than
62 or equal to 1, chosen so that the components of x will be less than
63 the overflow threshold. If the unscaled problem will not cause
64 overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A
65 is singular (A(j,j) = 0 for some j), then s is set to 0 and a
66 non-trivial solution to A*x = 0 is returned.
67
68 Arguments
69 =========
70
71 UPLO (input) CHARACTER*1
72 Specifies whether the matrix A is upper or lower triangular.
73 = 'U': Upper triangular
74 = 'L': Lower triangular
75
76 TRANS (input) CHARACTER*1
77 Specifies the operation applied to A.
78 = 'N': Solve A * x = s*b (No transpose)
79 = 'T': Solve A'* x = s*b (Transpose)
80 = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose)
81
82 DIAG (input) CHARACTER*1
83 Specifies whether or not the matrix A is unit triangular.
84 = 'N': Non-unit triangular
85 = 'U': Unit triangular
86
87 NORMIN (input) CHARACTER*1
88 Specifies whether CNORM has been set or not.
89 = 'Y': CNORM contains the column norms on entry
90 = 'N': CNORM is not set on entry. On exit, the norms will
91 be computed and stored in CNORM.
92
93 N (input) INTEGER
94 The order of the matrix A. N >= 0.
95
96 A (input) DOUBLE PRECISION array, dimension (LDA,N)
97 The triangular matrix A. If UPLO = 'U', the leading n by n
98 upper triangular part of the array A contains the upper
99 triangular matrix, and the strictly lower triangular part of
100 A is not referenced. If UPLO = 'L', the leading n by n lower
101 triangular part of the array A contains the lower triangular
102 matrix, and the strictly upper triangular part of A is not
103 referenced. If DIAG = 'U', the diagonal elements of A are
104 also not referenced and are assumed to be 1.
105
106 LDA (input) INTEGER
107 The leading dimension of the array A. LDA >= max (1,N).
108
109 X (input/output) DOUBLE PRECISION array, dimension (N)
110 On entry, the right hand side b of the triangular system.
111 On exit, X is overwritten by the solution vector x.
112
113 SCALE (output) DOUBLE PRECISION
114 The scaling factor s for the triangular system
115 A * x = s*b or A'* x = s*b.
116 If SCALE = 0, the matrix A is singular or badly scaled, and
117 the vector x is an exact or approximate solution to A*x = 0.
118
119 CNORM (input or output) DOUBLE PRECISION array, dimension (N)
120
121 If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
122 contains the norm of the off-diagonal part of the j-th column
123 of A. If TRANS = 'N', CNORM(j) must be greater than or equal
124 to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
125 must be greater than or equal to the 1-norm.
126
127 If NORMIN = 'N', CNORM is an output argument and CNORM(j)
128 returns the 1-norm of the offdiagonal part of the j-th column
129 of A.
130
131 INFO (output) INTEGER
132 = 0: successful exit
133 < 0: if INFO = -k, the k-th argument had an illegal value
134
135 Further Details
136 ======= =======
137
138 A rough bound on x is computed; if that is less than overflow, DTRSV
139 is called, otherwise, specific code is used which checks for possible
140 overflow or divide-by-zero at every operation.
141
142 A columnwise scheme is used for solving A*x = b. The basic algorithm
143 if A is lower triangular is
144
145 x[1:n] := b[1:n]
146 for j = 1, ..., n
147 x(j) := x(j) / A(j,j)
148 x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
149 end
150
151 Define bounds on the components of x after j iterations of the loop:
152 M(j) = bound on x[1:j]
153 G(j) = bound on x[j+1:n]
154 Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
155
156 Then for iteration j+1 we have
157 M(j+1) <= G(j) / | A(j+1,j+1) |
158 G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
159 <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
160
161 where CNORM(j+1) is greater than or equal to the infinity-norm of
162 column j+1 of A, not counting the diagonal. Hence
163
164 G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
165 1<=i<=j
166 and
167
168 |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
169 1<=i< j
170
171 Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the
172 reciprocal of the largest M(j), j=1,..,n, is larger than
173 max(underflow, 1/overflow).
174
175 The bound on x(j) is also used to determine when a step in the
176 columnwise method can be performed without fear of overflow. If
177 the computed bound is greater than a large constant, x is scaled to
178 prevent overflow, but if the bound overflows, x is set to 0, x(j) to
179 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
180
181 Similarly, a row-wise scheme is used to solve A'*x = b. The basic
182 algorithm for A upper triangular is
183
184 for j = 1, ..., n
185 x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
186 end
187
188 We simultaneously compute two bounds
189 G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
190 M(j) = bound on x(i), 1<=i<=j
191
192 The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
193 add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
194 Then the bound on x(j) is
195
196 M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
197
198 <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
199 1<=i<=j
200
201 and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater
202 than max(underflow, 1/overflow).
203
204 =====================================================================
205
206
207 Parameter adjustments */
208 /* Table of constant values */
209 integer c__1 = 1;
210 Treal c_b36 = .5;
211
212 /* System generated locals */
213 integer a_dim1, a_offset, i__1, i__2, i__3;
214 Treal d__1, d__2, d__3;
215 /* Local variables */
216 integer jinc;
217 Treal xbnd;
218 integer imax;
219 Treal tmax, tjjs, xmax, grow, sumj;
220 integer i__, j;
221 Treal tscal, uscal;
222 integer jlast;
223 logical upper;
224 Treal xj;
225 Treal bignum;
226 logical notran;
227 integer jfirst;
228 Treal smlnum;
229 logical nounit;
230 Treal rec, tjj;
231#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
232
233
234 a_dim1 = *lda;
235 a_offset = 1 + a_dim1 * 1;
236 a -= a_offset;
237 --x;
238 --cnorm;
239
240 /* Function Body */
241 *info = 0;
242 upper = template_blas_lsame(uplo, "U");
243 notran = template_blas_lsame(trans, "N");
244 nounit = template_blas_lsame(diag, "N");
245
246/* Test the input parameters. */
247
248 if (! upper && ! template_blas_lsame(uplo, "L")) {
249 *info = -1;
250 } else if (! notran && ! template_blas_lsame(trans, "T") && !
251 template_blas_lsame(trans, "C")) {
252 *info = -2;
253 } else if (! nounit && ! template_blas_lsame(diag, "U")) {
254 *info = -3;
255 } else if (! template_blas_lsame(normin, "Y") && ! template_blas_lsame(normin,
256 "N")) {
257 *info = -4;
258 } else if (*n < 0) {
259 *info = -5;
260 } else if (*lda < maxMACRO(1,*n)) {
261 *info = -7;
262 }
263 if (*info != 0) {
264 i__1 = -(*info);
265 template_blas_erbla("LATRS ", &i__1);
266 return 0;
267 }
268
269/* Quick return if possible */
270
271 if (*n == 0) {
272 return 0;
273 }
274
275/* Determine machine dependent parameters to control overflow. */
276
277 smlnum = template_lapack_lamch("Safe minimum", (Treal)0) / template_lapack_lamch("Precision", (Treal)0);
278 bignum = 1. / smlnum;
279 *scale = 1.;
280
281 if (template_blas_lsame(normin, "N")) {
282
283/* Compute the 1-norm of each column, not including the diagonal. */
284
285 if (upper) {
286
287/* A is upper triangular. */
288
289 i__1 = *n;
290 for (j = 1; j <= i__1; ++j) {
291 i__2 = j - 1;
292 cnorm[j] = template_blas_asum(&i__2, &a_ref(1, j), &c__1);
293/* L10: */
294 }
295 } else {
296
297/* A is lower triangular. */
298
299 i__1 = *n - 1;
300 for (j = 1; j <= i__1; ++j) {
301 i__2 = *n - j;
302 cnorm[j] = template_blas_asum(&i__2, &a_ref(j + 1, j), &c__1);
303/* L20: */
304 }
305 cnorm[*n] = 0.;
306 }
307 }
308
309/* Scale the column norms by TSCAL if the maximum element in CNORM is
310 greater than BIGNUM. */
311
312 imax = template_blas_idamax(n, &cnorm[1], &c__1);
313 tmax = cnorm[imax];
314 if (tmax <= bignum) {
315 tscal = 1.;
316 } else {
317 tscal = 1. / (smlnum * tmax);
318 dscal_(n, &tscal, &cnorm[1], &c__1);
319 }
320
321/* Compute a bound on the computed solution vector to see if the
322 Level 2 BLAS routine DTRSV can be used. */
323
324 j = template_blas_idamax(n, &x[1], &c__1);
325 xmax = (d__1 = x[j], absMACRO(d__1));
326 xbnd = xmax;
327 if (notran) {
328
329/* Compute the growth in A * x = b. */
330
331 if (upper) {
332 jfirst = *n;
333 jlast = 1;
334 jinc = -1;
335 } else {
336 jfirst = 1;
337 jlast = *n;
338 jinc = 1;
339 }
340
341 if (tscal != 1.) {
342 grow = 0.;
343 goto L50;
344 }
345
346 if (nounit) {
347
348/* A is non-unit triangular.
349
350 Compute GROW = 1/G(j) and XBND = 1/M(j).
351 Initially, G(0) = max{x(i), i=1,...,n}. */
352
353 grow = 1. / maxMACRO(xbnd,smlnum);
354 xbnd = grow;
355 i__1 = jlast;
356 i__2 = jinc;
357 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
358
359/* Exit the loop if the growth factor is too small. */
360
361 if (grow <= smlnum) {
362 goto L50;
363 }
364
365/* M(j) = G(j-1) / absMACRO(A(j,j)) */
366
367 tjj = (d__1 = a_ref(j, j), absMACRO(d__1));
368/* Computing MIN */
369 d__1 = xbnd, d__2 = minMACRO(1.,tjj) * grow;
370 xbnd = minMACRO(d__1,d__2);
371 if (tjj + cnorm[j] >= smlnum) {
372
373/* G(j) = G(j-1)*( 1 + CNORM(j) / absMACRO(A(j,j)) ) */
374
375 grow *= tjj / (tjj + cnorm[j]);
376 } else {
377
378/* G(j) could overflow, set GROW to 0. */
379
380 grow = 0.;
381 }
382/* L30: */
383 }
384 grow = xbnd;
385 } else {
386
387/* A is unit triangular.
388
389 Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
390
391 Computing MIN */
392 d__1 = 1., d__2 = 1. / maxMACRO(xbnd,smlnum);
393 grow = minMACRO(d__1,d__2);
394 i__2 = jlast;
395 i__1 = jinc;
396 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
397
398/* Exit the loop if the growth factor is too small. */
399
400 if (grow <= smlnum) {
401 goto L50;
402 }
403
404/* G(j) = G(j-1)*( 1 + CNORM(j) ) */
405
406 grow *= 1. / (cnorm[j] + 1.);
407/* L40: */
408 }
409 }
410L50:
411
412 ;
413 } else {
414
415/* Compute the growth in A' * x = b. */
416
417 if (upper) {
418 jfirst = 1;
419 jlast = *n;
420 jinc = 1;
421 } else {
422 jfirst = *n;
423 jlast = 1;
424 jinc = -1;
425 }
426
427 if (tscal != 1.) {
428 grow = 0.;
429 goto L80;
430 }
431
432 if (nounit) {
433
434/* A is non-unit triangular.
435
436 Compute GROW = 1/G(j) and XBND = 1/M(j).
437 Initially, M(0) = max{x(i), i=1,...,n}. */
438
439 grow = 1. / maxMACRO(xbnd,smlnum);
440 xbnd = grow;
441 i__1 = jlast;
442 i__2 = jinc;
443 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
444
445/* Exit the loop if the growth factor is too small. */
446
447 if (grow <= smlnum) {
448 goto L80;
449 }
450
451/* G(j) = maxMACRO( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) */
452
453 xj = cnorm[j] + 1.;
454/* Computing MIN */
455 d__1 = grow, d__2 = xbnd / xj;
456 grow = minMACRO(d__1,d__2);
457
458/* M(j) = M(j-1)*( 1 + CNORM(j) ) / absMACRO(A(j,j)) */
459
460 tjj = (d__1 = a_ref(j, j), absMACRO(d__1));
461 if (xj > tjj) {
462 xbnd *= tjj / xj;
463 }
464/* L60: */
465 }
466 grow = minMACRO(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 Computing MIN */
474 d__1 = 1., d__2 = 1. / maxMACRO(xbnd,smlnum);
475 grow = minMACRO(d__1,d__2);
476 i__2 = jlast;
477 i__1 = jinc;
478 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
479
480/* Exit the loop if the growth factor is too small. */
481
482 if (grow <= smlnum) {
483 goto L80;
484 }
485
486/* G(j) = ( 1 + CNORM(j) )*G(j-1) */
487
488 xj = cnorm[j] + 1.;
489 grow /= xj;
490/* L70: */
491 }
492 }
493L80:
494 ;
495 }
496
497 if (grow * tscal > smlnum) {
498
499/* Use the Level 2 BLAS solve if the reciprocal of the bound on
500 elements of X is not too small. */
501
502 template_blas_trsv(uplo, trans, diag, n, &a[a_offset], lda, &x[1], &c__1);
503 } else {
504
505/* Use a Level 1 BLAS solve, scaling intermediate results. */
506
507 if (xmax > bignum) {
508
509/* Scale X so that its components are less than or equal to
510 BIGNUM in absolute value. */
511
512 *scale = bignum / xmax;
513 dscal_(n, scale, &x[1], &c__1);
514 xmax = bignum;
515 }
516
517 if (notran) {
518
519/* Solve A * x = b */
520
521 i__1 = jlast;
522 i__2 = jinc;
523 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
524
525/* Compute x(j) = b(j) / A(j,j), scaling x if necessary. */
526
527 xj = (d__1 = x[j], absMACRO(d__1));
528 if (nounit) {
529 tjjs = a_ref(j, j) * tscal;
530 } else {
531 tjjs = tscal;
532 if (tscal == 1.) {
533 goto L100;
534 }
535 }
536 tjj = absMACRO(tjjs);
537 if (tjj > smlnum) {
538
539/* absMACRO(A(j,j)) > SMLNUM: */
540
541 if (tjj < 1.) {
542 if (xj > tjj * bignum) {
543
544/* Scale x by 1/b(j). */
545
546 rec = 1. / xj;
547 dscal_(n, &rec, &x[1], &c__1);
548 *scale *= rec;
549 xmax *= rec;
550 }
551 }
552 x[j] /= tjjs;
553 xj = (d__1 = x[j], absMACRO(d__1));
554 } else if (tjj > 0.) {
555
556/* 0 < absMACRO(A(j,j)) <= SMLNUM: */
557
558 if (xj > tjj * bignum) {
559
560/* Scale x by (1/absMACRO(x(j)))*absMACRO(A(j,j))*BIGNUM
561 to avoid overflow when dividing by A(j,j). */
562
563 rec = tjj * bignum / xj;
564 if (cnorm[j] > 1.) {
565
566/* Scale by 1/CNORM(j) to avoid overflow when
567 multiplying x(j) times column j. */
568
569 rec /= cnorm[j];
570 }
571 dscal_(n, &rec, &x[1], &c__1);
572 *scale *= rec;
573 xmax *= rec;
574 }
575 x[j] /= tjjs;
576 xj = (d__1 = x[j], absMACRO(d__1));
577 } else {
578
579/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
580 scale = 0, and compute a solution to A*x = 0. */
581
582 i__3 = *n;
583 for (i__ = 1; i__ <= i__3; ++i__) {
584 x[i__] = 0.;
585/* L90: */
586 }
587 x[j] = 1.;
588 xj = 1.;
589 *scale = 0.;
590 xmax = 0.;
591 }
592L100:
593
594/* Scale x if necessary to avoid overflow when adding a
595 multiple of column j of A. */
596
597 if (xj > 1.) {
598 rec = 1. / xj;
599 if (cnorm[j] > (bignum - xmax) * rec) {
600
601/* Scale x by 1/(2*absMACRO(x(j))). */
602
603 rec *= .5;
604 dscal_(n, &rec, &x[1], &c__1);
605 *scale *= rec;
606 }
607 } else if (xj * cnorm[j] > bignum - xmax) {
608
609/* Scale x by 1/2. */
610
611 dscal_(n, &c_b36, &x[1], &c__1);
612 *scale *= .5;
613 }
614
615 if (upper) {
616 if (j > 1) {
617
618/* Compute the update
619 x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) */
620
621 i__3 = j - 1;
622 d__1 = -x[j] * tscal;
623 daxpy_(&i__3, &d__1, &a_ref(1, j), &c__1, &x[1], &
624 c__1);
625 i__3 = j - 1;
626 i__ = template_blas_idamax(&i__3, &x[1], &c__1);
627 xmax = (d__1 = x[i__], absMACRO(d__1));
628 }
629 } else {
630 if (j < *n) {
631
632/* Compute the update
633 x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) */
634
635 i__3 = *n - j;
636 d__1 = -x[j] * tscal;
637 daxpy_(&i__3, &d__1, &a_ref(j + 1, j), &c__1, &x[j +
638 1], &c__1);
639 i__3 = *n - j;
640 i__ = j + template_blas_idamax(&i__3, &x[j + 1], &c__1);
641 xmax = (d__1 = x[i__], absMACRO(d__1));
642 }
643 }
644/* L110: */
645 }
646
647 } else {
648
649/* Solve A' * x = b */
650
651 i__2 = jlast;
652 i__1 = jinc;
653 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
654
655/* Compute x(j) = b(j) - sum A(k,j)*x(k).
656 k<>j */
657
658 xj = (d__1 = x[j], absMACRO(d__1));
659 uscal = tscal;
660 rec = 1. / maxMACRO(xmax,1.);
661 if (cnorm[j] > (bignum - xj) * rec) {
662
663/* If x(j) could overflow, scale x by 1/(2*XMAX). */
664
665 rec *= .5;
666 if (nounit) {
667 tjjs = a_ref(j, j) * tscal;
668 } else {
669 tjjs = tscal;
670 }
671 tjj = absMACRO(tjjs);
672 if (tjj > 1.) {
673
674/* Divide by A(j,j) when scaling x if A(j,j) > 1.
675
676 Computing MIN */
677 d__1 = 1., d__2 = rec * tjj;
678 rec = minMACRO(d__1,d__2);
679 uscal /= tjjs;
680 }
681 if (rec < 1.) {
682 dscal_(n, &rec, &x[1], &c__1);
683 *scale *= rec;
684 xmax *= rec;
685 }
686 }
687
688 sumj = 0.;
689 if (uscal == 1.) {
690
691/* If the scaling needed for A in the dot product is 1,
692 call DDOT to perform the dot product. */
693
694 if (upper) {
695 i__3 = j - 1;
696 sumj = ddot_(&i__3, &a_ref(1, j), &c__1, &x[1], &c__1)
697 ;
698 } else if (j < *n) {
699 i__3 = *n - j;
700 sumj = ddot_(&i__3, &a_ref(j + 1, j), &c__1, &x[j + 1]
701 , &c__1);
702 }
703 } else {
704
705/* Otherwise, use in-line code for the dot product. */
706
707 if (upper) {
708 i__3 = j - 1;
709 for (i__ = 1; i__ <= i__3; ++i__) {
710 sumj += a_ref(i__, j) * uscal * x[i__];
711/* L120: */
712 }
713 } else if (j < *n) {
714 i__3 = *n;
715 for (i__ = j + 1; i__ <= i__3; ++i__) {
716 sumj += a_ref(i__, j) * uscal * x[i__];
717/* L130: */
718 }
719 }
720 }
721
722 if (uscal == tscal) {
723
724/* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
725 was not used to scale the dotproduct. */
726
727 x[j] -= sumj;
728 xj = (d__1 = x[j], absMACRO(d__1));
729 if (nounit) {
730 tjjs = a_ref(j, j) * tscal;
731 } else {
732 tjjs = tscal;
733 if (tscal == 1.) {
734 goto L150;
735 }
736 }
737
738/* Compute x(j) = x(j) / A(j,j), scaling if necessary. */
739
740 tjj = absMACRO(tjjs);
741 if (tjj > smlnum) {
742
743/* absMACRO(A(j,j)) > SMLNUM: */
744
745 if (tjj < 1.) {
746 if (xj > tjj * bignum) {
747
748/* Scale X by 1/absMACRO(x(j)). */
749
750 rec = 1. / xj;
751 dscal_(n, &rec, &x[1], &c__1);
752 *scale *= rec;
753 xmax *= rec;
754 }
755 }
756 x[j] /= tjjs;
757 } else if (tjj > 0.) {
758
759/* 0 < absMACRO(A(j,j)) <= SMLNUM: */
760
761 if (xj > tjj * bignum) {
762
763/* Scale x by (1/absMACRO(x(j)))*absMACRO(A(j,j))*BIGNUM. */
764
765 rec = tjj * bignum / xj;
766 dscal_(n, &rec, &x[1], &c__1);
767 *scale *= rec;
768 xmax *= rec;
769 }
770 x[j] /= tjjs;
771 } else {
772
773/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
774 scale = 0, and compute a solution to A'*x = 0. */
775
776 i__3 = *n;
777 for (i__ = 1; i__ <= i__3; ++i__) {
778 x[i__] = 0.;
779/* L140: */
780 }
781 x[j] = 1.;
782 *scale = 0.;
783 xmax = 0.;
784 }
785L150:
786 ;
787 } else {
788
789/* Compute x(j) := x(j) / A(j,j) - sumj if the dot
790 product has already been divided by 1/A(j,j). */
791
792 x[j] = x[j] / tjjs - sumj;
793 }
794/* Computing MAX */
795 d__2 = xmax, d__3 = (d__1 = x[j], absMACRO(d__1));
796 xmax = maxMACRO(d__2,d__3);
797/* L160: */
798 }
799 }
800 *scale /= tscal;
801 }
802
803/* Scale the column norms by 1/TSCAL for return. */
804
805 if (tscal != 1.) {
806 d__1 = 1. / tscal;
807 dscal_(n, &d__1, &cnorm[1], &c__1);
808 }
809
810 return 0;
811
812/* End of DLATRS */
813
814} /* dlatrs_ */
815
816#undef a_ref
817
818
819#endif
void dscal_(const int *n, const double *da, double *dx, const int *incx)
double ddot_(const int *n, const double *dx, const int *incx, const double *dy, const int *incy)
void daxpy_(const int *n, const double *da, const double *dx, const int *incx, double *dy, const int *incy)
Treal template_blas_asum(const integer *n, const Treal *dx, const integer *incx)
Definition template_blas_asum.h:42
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition template_blas_common.cc:46
int integer
Definition template_blas_common.h:40
#define absMACRO(x)
Definition template_blas_common.h:47
#define minMACRO(a, b)
Definition template_blas_common.h:46
#define maxMACRO(a, b)
Definition template_blas_common.h:45
bool logical
Definition template_blas_common.h:41
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition template_blas_idamax.h:42
int template_blas_trsv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *a, const integer *lda, Treal *x, const integer *incx)
Definition template_blas_trsv.h:42
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition template_lapack_lamch.h:202
int template_lapack_latrs(const char *uplo, const char *trans, const char *diag, const char *normin, const integer *n, const Treal *a, const integer *lda, Treal *x, Treal *scale, Treal *cnorm, integer *info)
Definition template_lapack_latrs.h:42
#define a_ref(a_1, a_2)