ergo
template_lapack_larfg.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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_LARFG_HEADER
38 #define TEMPLATE_LAPACK_LARFG_HEADER
39 
40 #include "template_lapack_lapy2.h"
41 
42 template<class Treal>
43 int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x,
44  const integer *incx, Treal *tau)
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  September 30, 1994
50 
51 
52  Purpose
53  =======
54 
55  DLARFG generates a real elementary reflector H of order n, such
56  that
57 
58  H * ( alpha ) = ( beta ), H' * H = I.
59  ( x ) ( 0 )
60 
61  where alpha and beta are scalars, and x is an (n-1)-element real
62  vector. H is represented in the form
63 
64  H = I - tau * ( 1 ) * ( 1 v' ) ,
65  ( v )
66 
67  where tau is a real scalar and v is a real (n-1)-element
68  vector.
69 
70  If the elements of x are all zero, then tau = 0 and H is taken to be
71  the unit matrix.
72 
73  Otherwise 1 <= tau <= 2.
74 
75  Arguments
76  =========
77 
78  N (input) INTEGER
79  The order of the elementary reflector.
80 
81  ALPHA (input/output) DOUBLE PRECISION
82  On entry, the value alpha.
83  On exit, it is overwritten with the value beta.
84 
85  X (input/output) DOUBLE PRECISION array, dimension
86  (1+(N-2)*abs(INCX))
87  On entry, the vector x.
88  On exit, it is overwritten with the vector v.
89 
90  INCX (input) INTEGER
91  The increment between elements of X. INCX > 0.
92 
93  TAU (output) DOUBLE PRECISION
94  The value tau.
95 
96  =====================================================================
97 
98 
99  Parameter adjustments */
100  /* System generated locals */
101  integer i__1;
102  Treal d__1;
103  /* Local variables */
104  Treal beta;
105  integer j;
106  Treal xnorm;
107  Treal safmin, rsafmn;
108  integer knt;
109 
110  --x;
111 
112  /* Function Body */
113  if (*n <= 1) {
114  *tau = 0.;
115  return 0;
116  }
117 
118  i__1 = *n - 1;
119  xnorm = template_blas_nrm2(&i__1, &x[1], incx);
120 
121  if (xnorm == 0.) {
122 
123 /* H = I */
124 
125  *tau = 0.;
126  } else {
127 
128 /* general case */
129 
130  d__1 = template_lapack_lapy2(alpha, &xnorm);
131  beta = -template_lapack_d_sign(&d__1, alpha);
132  safmin = template_lapack_lamch("S", (Treal)0) / template_lapack_lamch("E", (Treal)0);
133  if (absMACRO(beta) < safmin) {
134 
135 /* XNORM, BETA may be inaccurate; scale X and recompute them */
136 
137  rsafmn = 1. / safmin;
138  knt = 0;
139 L10:
140  ++knt;
141  i__1 = *n - 1;
142  template_blas_scal(&i__1, &rsafmn, &x[1], incx);
143  beta *= rsafmn;
144  *alpha *= rsafmn;
145  if (absMACRO(beta) < safmin) {
146  goto L10;
147  }
148 
149 /* New BETA is at most 1, at least SAFMIN */
150 
151  i__1 = *n - 1;
152  xnorm = template_blas_nrm2(&i__1, &x[1], incx);
153  d__1 = template_lapack_lapy2(alpha, &xnorm);
154  beta = -template_lapack_d_sign(&d__1, alpha);
155  *tau = (beta - *alpha) / beta;
156  i__1 = *n - 1;
157  d__1 = 1. / (*alpha - beta);
158  template_blas_scal(&i__1, &d__1, &x[1], incx);
159 
160 /* If ALPHA is subnormal, it may lose relative accuracy */
161 
162  *alpha = beta;
163  i__1 = knt;
164  for (j = 1; j <= i__1; ++j) {
165  *alpha *= safmin;
166 /* L20: */
167  }
168  } else {
169  *tau = (beta - *alpha) / beta;
170  i__1 = *n - 1;
171  d__1 = 1. / (*alpha - beta);
172  template_blas_scal(&i__1, &d__1, &x[1], incx);
173  *alpha = beta;
174  }
175  }
176 
177  return 0;
178 
179 /* End of DLARFG */
180 
181 } /* dlarfg_ */
182 
183 #endif
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:43
Treal template_blas_nrm2(const integer *n, const Treal *x, const integer *incx)
Definition: template_blas_nrm2.h:42
#define absMACRO(x)
Definition: template_blas_common.h:47
int integer
Definition: template_blas_common.h:40
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:42
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:48
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition: template_lapack_larfg.h:43
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202