00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_LASV2_HEADER 00038 #define TEMPLATE_LAPACK_LASV2_HEADER 00039 00040 00041 template<class Treal> 00042 int template_lapack_lasv2(const Treal *f, const Treal *g, const Treal *h__, 00043 Treal *ssmin, Treal *ssmax, Treal *snr, Treal * 00044 csr, Treal *snl, Treal *csl) 00045 { 00046 /* -- LAPACK auxiliary routine (version 3.0) -- 00047 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00048 Courant Institute, Argonne National Lab, and Rice University 00049 October 31, 1992 00050 00051 00052 Purpose 00053 ======= 00054 00055 DLASV2 computes the singular value decomposition of a 2-by-2 00056 triangular matrix 00057 [ F G ] 00058 [ 0 H ]. 00059 On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the 00060 smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and 00061 right singular vectors for abs(SSMAX), giving the decomposition 00062 00063 [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ] 00064 [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ]. 00065 00066 Arguments 00067 ========= 00068 00069 F (input) DOUBLE PRECISION 00070 The (1,1) element of the 2-by-2 matrix. 00071 00072 G (input) DOUBLE PRECISION 00073 The (1,2) element of the 2-by-2 matrix. 00074 00075 H (input) DOUBLE PRECISION 00076 The (2,2) element of the 2-by-2 matrix. 00077 00078 SSMIN (output) DOUBLE PRECISION 00079 abs(SSMIN) is the smaller singular value. 00080 00081 SSMAX (output) DOUBLE PRECISION 00082 abs(SSMAX) is the larger singular value. 00083 00084 SNL (output) DOUBLE PRECISION 00085 CSL (output) DOUBLE PRECISION 00086 The vector (CSL, SNL) is a unit left singular vector for the 00087 singular value abs(SSMAX). 00088 00089 SNR (output) DOUBLE PRECISION 00090 CSR (output) DOUBLE PRECISION 00091 The vector (CSR, SNR) is a unit right singular vector for the 00092 singular value abs(SSMAX). 00093 00094 Further Details 00095 =============== 00096 00097 Any input parameter may be aliased with any output parameter. 00098 00099 Barring over/underflow and assuming a guard digit in subtraction, all 00100 output quantities are correct to within a few units in the last 00101 place (ulps). 00102 00103 In IEEE arithmetic, the code works correctly if one matrix element is 00104 infinite. 00105 00106 Overflow will not occur unless the largest singular value itself 00107 overflows or is within a few ulps of overflow. (On machines with 00108 partial overflow, like the Cray, overflow may occur if the largest 00109 singular value is within a factor of 2 of overflow.) 00110 00111 Underflow is harmless if underflow is gradual. Otherwise, results 00112 may correspond to a matrix modified by perturbations of size near 00113 the underflow threshold. 00114 00115 ===================================================================== */ 00116 /* Table of constant values */ 00117 Treal c_b3 = 2.; 00118 Treal c_b4 = 1.; 00119 00120 /* System generated locals */ 00121 Treal d__1; 00122 /* Local variables */ 00123 integer pmax; 00124 Treal temp; 00125 logical swap; 00126 Treal a, d__, l, m, r__, s, t, tsign, fa, ga, ha; 00127 Treal ft, gt, ht, mm; 00128 logical gasmal; 00129 Treal tt, clt, crt, slt, srt; 00130 00131 /* Initialization added by Elias to get rid of compiler warnings. */ 00132 tsign = 0; 00133 00134 00135 ft = *f; 00136 fa = absMACRO(ft); 00137 ht = *h__; 00138 ha = absMACRO(*h__); 00139 00140 /* PMAX points to the maximum absolute element of matrix 00141 PMAX = 1 if F largest in absolute values 00142 PMAX = 2 if G largest in absolute values 00143 PMAX = 3 if H largest in absolute values */ 00144 00145 pmax = 1; 00146 swap = ha > fa; 00147 if (swap) { 00148 pmax = 3; 00149 temp = ft; 00150 ft = ht; 00151 ht = temp; 00152 temp = fa; 00153 fa = ha; 00154 ha = temp; 00155 00156 /* Now FA .ge. HA */ 00157 00158 } 00159 gt = *g; 00160 ga = absMACRO(gt); 00161 if (ga == 0.) { 00162 00163 /* Diagonal matrix */ 00164 00165 *ssmin = ha; 00166 *ssmax = fa; 00167 clt = 1.; 00168 crt = 1.; 00169 slt = 0.; 00170 srt = 0.; 00171 } else { 00172 gasmal = TRUE_; 00173 if (ga > fa) { 00174 pmax = 2; 00175 if (fa / ga < template_lapack_lamch("EPS", (Treal)0)) { 00176 00177 /* Case of very large GA */ 00178 00179 gasmal = FALSE_; 00180 *ssmax = ga; 00181 if (ha > 1.) { 00182 *ssmin = fa / (ga / ha); 00183 } else { 00184 *ssmin = fa / ga * ha; 00185 } 00186 clt = 1.; 00187 slt = ht / gt; 00188 srt = 1.; 00189 crt = ft / gt; 00190 } 00191 } 00192 if (gasmal) { 00193 00194 /* Normal case */ 00195 00196 d__ = fa - ha; 00197 if (d__ == fa) { 00198 00199 /* Copes with infinite F or H */ 00200 00201 l = 1.; 00202 } else { 00203 l = d__ / fa; 00204 } 00205 00206 /* Note that 0 .le. L .le. 1 */ 00207 00208 m = gt / ft; 00209 00210 /* Note that abs(M) .le. 1/macheps */ 00211 00212 t = 2. - l; 00213 00214 /* Note that T .ge. 1 */ 00215 00216 mm = m * m; 00217 tt = t * t; 00218 s = template_blas_sqrt(tt + mm); 00219 00220 /* Note that 1 .le. S .le. 1 + 1/macheps */ 00221 00222 if (l == 0.) { 00223 r__ = absMACRO(m); 00224 } else { 00225 r__ = template_blas_sqrt(l * l + mm); 00226 } 00227 00228 /* Note that 0 .le. R .le. 1 + 1/macheps */ 00229 00230 a = (s + r__) * .5; 00231 00232 /* Note that 1 .le. A .le. 1 + abs(M) */ 00233 00234 *ssmin = ha / a; 00235 *ssmax = fa * a; 00236 if (mm == 0.) { 00237 00238 /* Note that M is very tiny */ 00239 00240 if (l == 0.) { 00241 t = template_lapack_d_sign(&c_b3, &ft) * template_lapack_d_sign(&c_b4, >); 00242 } else { 00243 t = gt / template_lapack_d_sign(&d__, &ft) + m / t; 00244 } 00245 } else { 00246 t = (m / (s + t) + m / (r__ + l)) * (a + 1.); 00247 } 00248 l = template_blas_sqrt(t * t + 4.); 00249 crt = 2. / l; 00250 srt = t / l; 00251 clt = (crt + srt * m) / a; 00252 slt = ht / ft * srt / a; 00253 } 00254 } 00255 if (swap) { 00256 *csl = srt; 00257 *snl = crt; 00258 *csr = slt; 00259 *snr = clt; 00260 } else { 00261 *csl = clt; 00262 *snl = slt; 00263 *csr = crt; 00264 *snr = srt; 00265 } 00266 00267 /* Correct signs of SSMAX and SSMIN */ 00268 00269 if (pmax == 1) { 00270 tsign = template_lapack_d_sign(&c_b4, csr) * template_lapack_d_sign(&c_b4, csl) * template_lapack_d_sign(&c_b4, f); 00271 } 00272 if (pmax == 2) { 00273 tsign = template_lapack_d_sign(&c_b4, snr) * template_lapack_d_sign(&c_b4, csl) * template_lapack_d_sign(&c_b4, g); 00274 } 00275 if (pmax == 3) { 00276 tsign = template_lapack_d_sign(&c_b4, snr) * template_lapack_d_sign(&c_b4, snl) * template_lapack_d_sign(&c_b4, h__); 00277 } 00278 *ssmax = template_lapack_d_sign(ssmax, &tsign); 00279 d__1 = tsign * template_lapack_d_sign(&c_b4, f) * template_lapack_d_sign(&c_b4, h__); 00280 *ssmin = template_lapack_d_sign(ssmin, &d__1); 00281 return 0; 00282 00283 /* End of DLASV2 */ 00284 00285 } /* dlasv2_ */ 00286 00287 #endif