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