ergo
template_lapack_lartg.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_LARTG_HEADER
38#define TEMPLATE_LAPACK_LARTG_HEADER
39
40
41template<class Treal>
42int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs,
43 Treal *sn, Treal *r__)
44{
45/* -- LAPACK auxiliary routine (version 3.0) --
46 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
47 Courant Institute, Argonne National Lab, and Rice University
48 September 30, 1994
49
50
51 Purpose
52 =======
53
54 DLARTG generate a plane rotation so that
55
56 [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1.
57 [ -SN CS ] [ G ] [ 0 ]
58
59 This is a slower, more accurate version of the BLAS1 routine DROTG,
60 with the following other differences:
61 F and G are unchanged on return.
62 If G=0, then CS=1 and SN=0.
63 If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any
64 floating point operations (saves work in DBDSQR when
65 there are zeros on the diagonal).
66
67 If F exceeds G in magnitude, CS will be positive.
68
69 Arguments
70 =========
71
72 F (input) DOUBLE PRECISION
73 The first component of vector to be rotated.
74
75 G (input) DOUBLE PRECISION
76 The second component of vector to be rotated.
77
78 CS (output) DOUBLE PRECISION
79 The cosine of the rotation.
80
81 SN (output) DOUBLE PRECISION
82 The sine of the rotation.
83
84 R (output) DOUBLE PRECISION
85 The nonzero component of the rotated vector.
86
87 ===================================================================== */
88 /* Initialized data */
89 logical first = TRUE_;
90 /* System generated locals */
91 integer i__1;
92 Treal d__1, d__2;
93 /* Local variables */
94 integer i__;
95 Treal scale;
96 integer count;
97 Treal f1, g1, safmn2, safmx2;
98 Treal safmin, eps;
99
100
101
102 if (first) {
103 first = FALSE_;
104 safmin = template_lapack_lamch("S", (Treal)0);
105 eps = template_lapack_lamch("E", (Treal)0);
106 d__1 = template_lapack_lamch("B", (Treal)0);
107 i__1 = (integer) (template_blas_log(safmin / eps) / template_blas_log(template_lapack_lamch("B", (Treal)0)) /
108 2.);
109 safmn2 = template_lapack_pow_di(&d__1, &i__1);
110 safmx2 = 1. / safmn2;
111 }
112 if (*g == 0.) {
113 *cs = 1.;
114 *sn = 0.;
115 *r__ = *f;
116 } else if (*f == 0.) {
117 *cs = 0.;
118 *sn = 1.;
119 *r__ = *g;
120 } else {
121 f1 = *f;
122 g1 = *g;
123/* Computing MAX */
124 d__1 = absMACRO(f1), d__2 = absMACRO(g1);
125 scale = maxMACRO(d__1,d__2);
126 if (scale >= safmx2) {
127 count = 0;
128L10:
129 ++count;
130 f1 *= safmn2;
131 g1 *= safmn2;
132/* Computing MAX */
133 d__1 = absMACRO(f1), d__2 = absMACRO(g1);
134 scale = maxMACRO(d__1,d__2);
135 if (scale >= safmx2) {
136 goto L10;
137 }
138/* Computing 2nd power */
139 d__1 = f1;
140/* Computing 2nd power */
141 d__2 = g1;
142 *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2);
143 *cs = f1 / *r__;
144 *sn = g1 / *r__;
145 i__1 = count;
146 for (i__ = 1; i__ <= i__1; ++i__) {
147 *r__ *= safmx2;
148/* L20: */
149 }
150 } else if (scale <= safmn2) {
151 count = 0;
152L30:
153 ++count;
154 f1 *= safmx2;
155 g1 *= safmx2;
156/* Computing MAX */
157 d__1 = absMACRO(f1), d__2 = absMACRO(g1);
158 scale = maxMACRO(d__1,d__2);
159 if (scale <= safmn2) {
160 goto L30;
161 }
162/* Computing 2nd power */
163 d__1 = f1;
164/* Computing 2nd power */
165 d__2 = g1;
166 *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2);
167 *cs = f1 / *r__;
168 *sn = g1 / *r__;
169 i__1 = count;
170 for (i__ = 1; i__ <= i__1; ++i__) {
171 *r__ *= safmn2;
172/* L40: */
173 }
174 } else {
175/* Computing 2nd power */
176 d__1 = f1;
177/* Computing 2nd power */
178 d__2 = g1;
179 *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2);
180 *cs = f1 / *r__;
181 *sn = g1 / *r__;
182 }
183 if (absMACRO(*f) > absMACRO(*g) && *cs < 0.) {
184 *cs = -(*cs);
185 *sn = -(*sn);
186 *r__ = -(*r__);
187 }
188 }
189 return 0;
190
191/* End of DLARTG */
192
193} /* dlartg_ */
194
195#endif
Treal template_blas_sqrt(Treal x)
Treal template_blas_log(Treal x)
int integer
Definition: template_blas_common.h:40
#define absMACRO(x)
Definition: template_blas_common.h:47
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
#define TRUE_
Definition: template_lapack_common.h:42
#define FALSE_
Definition: template_lapack_common.h:43
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
double template_lapack_pow_di(Treal *ap, integer *bp)
Definition: template_lapack_lamch.h:168
int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition: template_lapack_lartg.h:42