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