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