ergo
template_lapack_lartg.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_LARTG_HEADER
38 #define TEMPLATE_LAPACK_LARTG_HEADER
39 
40 
41 template<class Treal>
42 int 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;
128 L10:
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;
152 L30:
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
#define absMACRO(x)
Definition: template_blas_common.h:47
int integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
Treal template_blas_log(Treal x)
int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition: template_lapack_lartg.h:42
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
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_blas_sqrt(Treal x)