ergo
template_lapack_larfg.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_LARFG_HEADER
38#define TEMPLATE_LAPACK_LARFG_HEADER
39
41
42template<class Treal>
43int 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;
139L10:
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 integer
Definition: template_blas_common.h:40
#define absMACRO(x)
Definition: template_blas_common.h:47
Treal template_blas_nrm2(const integer *n, const Treal *x, const integer *incx)
Definition: template_blas_nrm2.h:42
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:43
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:48
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:42
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition: template_lapack_larfg.h:43