ergo
template_lapack_lasv2.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_LASV2_HEADER
38 #define TEMPLATE_LAPACK_LASV2_HEADER
39 
40 
41 template<class Treal>
42 int 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
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
template_lapack_lamch
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
absMACRO
#define absMACRO(x)
Definition: template_blas_common.h:47
logical
bool logical
Definition: template_blas_common.h:41
TRUE_
#define TRUE_
Definition: template_lapack_common.h:42
integer
int integer
Definition: template_blas_common.h:40
template_lapack_d_sign
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:48
FALSE_
#define FALSE_
Definition: template_lapack_common.h:43
template_lapack_lasv2
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