ergo
template_lapack_lasv2.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_LASV2_HEADER
38#define TEMPLATE_LAPACK_LASV2_HEADER
39
40
41template<class Treal>
42int template_lapack_lasv2(const Treal *f, const Treal *g, const Treal *h__,
43 Treal *ssmin, Treal *ssmax, Treal *snr, Treal *
44 csr, Treal *snl, Treal *csl)
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 October 31, 1992
50
51
52 Purpose
53 =======
54
55 DLASV2 computes the singular value decomposition of a 2-by-2
56 triangular matrix
57 [ F G ]
58 [ 0 H ].
59 On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
60 smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
61 right singular vectors for abs(SSMAX), giving the decomposition
62
63 [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ]
64 [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ].
65
66 Arguments
67 =========
68
69 F (input) DOUBLE PRECISION
70 The (1,1) element of the 2-by-2 matrix.
71
72 G (input) DOUBLE PRECISION
73 The (1,2) element of the 2-by-2 matrix.
74
75 H (input) DOUBLE PRECISION
76 The (2,2) element of the 2-by-2 matrix.
77
78 SSMIN (output) DOUBLE PRECISION
79 abs(SSMIN) is the smaller singular value.
80
81 SSMAX (output) DOUBLE PRECISION
82 abs(SSMAX) is the larger singular value.
83
84 SNL (output) DOUBLE PRECISION
85 CSL (output) DOUBLE PRECISION
86 The vector (CSL, SNL) is a unit left singular vector for the
87 singular value abs(SSMAX).
88
89 SNR (output) DOUBLE PRECISION
90 CSR (output) DOUBLE PRECISION
91 The vector (CSR, SNR) is a unit right singular vector for the
92 singular value abs(SSMAX).
93
94 Further Details
95 ===============
96
97 Any input parameter may be aliased with any output parameter.
98
99 Barring over/underflow and assuming a guard digit in subtraction, all
100 output quantities are correct to within a few units in the last
101 place (ulps).
102
103 In IEEE arithmetic, the code works correctly if one matrix element is
104 infinite.
105
106 Overflow will not occur unless the largest singular value itself
107 overflows or is within a few ulps of overflow. (On machines with
108 partial overflow, like the Cray, overflow may occur if the largest
109 singular value is within a factor of 2 of overflow.)
110
111 Underflow is harmless if underflow is gradual. Otherwise, results
112 may correspond to a matrix modified by perturbations of size near
113 the underflow threshold.
114
115 ===================================================================== */
116 /* Table of constant values */
117 Treal c_b3 = 2.;
118 Treal c_b4 = 1.;
119
120 /* System generated locals */
121 Treal d__1;
122 /* Local variables */
123 integer pmax;
124 Treal temp;
125 logical swap;
126 Treal a, d__, l, m, r__, s, t, tsign, fa, ga, ha;
127 Treal ft, gt, ht, mm;
128 logical gasmal;
129 Treal tt, clt, crt, slt, srt;
130
131 /* Initialization added by Elias to get rid of compiler warnings. */
132 tsign = 0;
133
134
135 ft = *f;
136 fa = absMACRO(ft);
137 ht = *h__;
138 ha = absMACRO(*h__);
139
140/* PMAX points to the maximum absolute element of matrix
141 PMAX = 1 if F largest in absolute values
142 PMAX = 2 if G largest in absolute values
143 PMAX = 3 if H largest in absolute values */
144
145 pmax = 1;
146 swap = ha > fa;
147 if (swap) {
148 pmax = 3;
149 temp = ft;
150 ft = ht;
151 ht = temp;
152 temp = fa;
153 fa = ha;
154 ha = temp;
155
156/* Now FA .ge. HA */
157
158 }
159 gt = *g;
160 ga = absMACRO(gt);
161 if (ga == 0.) {
162
163/* Diagonal matrix */
164
165 *ssmin = ha;
166 *ssmax = fa;
167 clt = 1.;
168 crt = 1.;
169 slt = 0.;
170 srt = 0.;
171 } else {
172 gasmal = TRUE_;
173 if (ga > fa) {
174 pmax = 2;
175 if (fa / ga < template_lapack_lamch("EPS", (Treal)0)) {
176
177/* Case of very large GA */
178
179 gasmal = FALSE_;
180 *ssmax = ga;
181 if (ha > 1.) {
182 *ssmin = fa / (ga / ha);
183 } else {
184 *ssmin = fa / ga * ha;
185 }
186 clt = 1.;
187 slt = ht / gt;
188 srt = 1.;
189 crt = ft / gt;
190 }
191 }
192 if (gasmal) {
193
194/* Normal case */
195
196 d__ = fa - ha;
197 if (d__ == fa) {
198
199/* Copes with infinite F or H */
200
201 l = 1.;
202 } else {
203 l = d__ / fa;
204 }
205
206/* Note that 0 .le. L .le. 1 */
207
208 m = gt / ft;
209
210/* Note that abs(M) .le. 1/macheps */
211
212 t = 2. - l;
213
214/* Note that T .ge. 1 */
215
216 mm = m * m;
217 tt = t * t;
218 s = template_blas_sqrt(tt + mm);
219
220/* Note that 1 .le. S .le. 1 + 1/macheps */
221
222 if (l == 0.) {
223 r__ = absMACRO(m);
224 } else {
225 r__ = template_blas_sqrt(l * l + mm);
226 }
227
228/* Note that 0 .le. R .le. 1 + 1/macheps */
229
230 a = (s + r__) * .5;
231
232/* Note that 1 .le. A .le. 1 + abs(M) */
233
234 *ssmin = ha / a;
235 *ssmax = fa * a;
236 if (mm == 0.) {
237
238/* Note that M is very tiny */
239
240 if (l == 0.) {
241 t = template_lapack_d_sign(&c_b3, &ft) * template_lapack_d_sign(&c_b4, &gt);
242 } else {
243 t = gt / template_lapack_d_sign(&d__, &ft) + m / t;
244 }
245 } else {
246 t = (m / (s + t) + m / (r__ + l)) * (a + 1.);
247 }
248 l = template_blas_sqrt(t * t + 4.);
249 crt = 2. / l;
250 srt = t / l;
251 clt = (crt + srt * m) / a;
252 slt = ht / ft * srt / a;
253 }
254 }
255 if (swap) {
256 *csl = srt;
257 *snl = crt;
258 *csr = slt;
259 *snr = clt;
260 } else {
261 *csl = clt;
262 *snl = slt;
263 *csr = crt;
264 *snr = srt;
265 }
266
267/* Correct signs of SSMAX and SSMIN */
268
269 if (pmax == 1) {
270 tsign = template_lapack_d_sign(&c_b4, csr) * template_lapack_d_sign(&c_b4, csl) * template_lapack_d_sign(&c_b4, f);
271 }
272 if (pmax == 2) {
273 tsign = template_lapack_d_sign(&c_b4, snr) * template_lapack_d_sign(&c_b4, csl) * template_lapack_d_sign(&c_b4, g);
274 }
275 if (pmax == 3) {
276 tsign = template_lapack_d_sign(&c_b4, snr) * template_lapack_d_sign(&c_b4, snl) * template_lapack_d_sign(&c_b4, h__);
277 }
278 *ssmax = template_lapack_d_sign(ssmax, &tsign);
279 d__1 = tsign * template_lapack_d_sign(&c_b4, f) * template_lapack_d_sign(&c_b4, h__);
280 *ssmin = template_lapack_d_sign(ssmin, &d__1);
281 return 0;
282
283/* End of DLASV2 */
284
285} /* dlasv2_ */
286
287#endif
Treal template_blas_sqrt(Treal x)
int integer
Definition: template_blas_common.h:40
#define absMACRO(x)
Definition: template_blas_common.h:47
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
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:48
int template_lapack_lasv2(const Treal *f, const Treal *g, const Treal *h__, Treal *ssmin, Treal *ssmax, Treal *snr, Treal *csr, Treal *snl, Treal *csl)
Definition: template_lapack_lasv2.h:42