ergo
template_lapack_latrs.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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 
41 template<class Treal>
42 int 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  }
410 L50:
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  }
493 L80:
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  }
592 L100:
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  }
785 L150:
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
template_blas_idamax
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_idamax.h:42
template_lapack_lamch
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
minMACRO
#define minMACRO(a, b)
Definition: template_blas_common.h:46
absMACRO
#define absMACRO(x)
Definition: template_blas_common.h:47
dscal_
void dscal_(const int *n, const double *da, double *dx, const int *incx)
ddot_
double ddot_(const int *n, const double *dx, const int *incx, const double *dy, const int *incy)
logical
bool logical
Definition: template_blas_common.h:41
template_blas_asum
Treal template_blas_asum(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_asum.h:42
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
integer
int integer
Definition: template_blas_common.h:40
template_lapack_latrs
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
a_ref
#define a_ref(a_1, a_2)
template_blas_trsv
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
daxpy_
void daxpy_(const int *n, const double *da, const double *dx, const int *incx, double *dy, const int *incy)
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45