ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00028 /* This file belongs to the template_lapack part of the Ergo source 00029 * code. The source files in the template_lapack directory are modified 00030 * versions of files originally distributed as CLAPACK, see the 00031 * Copyright/license notice in the file template_lapack/COPYING. 00032 */ 00033 00034 00035 #ifndef TEMPLATE_LAPACK_STERF_HEADER 00036 #define TEMPLATE_LAPACK_STERF_HEADER 00037 00038 #include "template_lapack_common.h" 00039 00040 template<class Treal> 00041 int template_lapack_sterf(const integer *n, Treal *d__, Treal *e, 00042 integer *info) 00043 { 00044 /* -- LAPACK routine (version 3.0) -- 00045 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00046 Courant Institute, Argonne National Lab, and Rice University 00047 June 30, 1999 00048 00049 00050 Purpose 00051 ======= 00052 00053 DSTERF computes all eigenvalues of a symmetric tridiagonal matrix 00054 using the Pal-Walker-Kahan variant of the QL or QR algorithm. 00055 00056 Arguments 00057 ========= 00058 00059 N (input) INTEGER 00060 The order of the matrix. N >= 0. 00061 00062 D (input/output) DOUBLE PRECISION array, dimension (N) 00063 On entry, the n diagonal elements of the tridiagonal matrix. 00064 On exit, if INFO = 0, the eigenvalues in ascending order. 00065 00066 E (input/output) DOUBLE PRECISION array, dimension (N-1) 00067 On entry, the (n-1) subdiagonal elements of the tridiagonal 00068 matrix. 00069 On exit, E has been destroyed. 00070 00071 INFO (output) INTEGER 00072 = 0: successful exit 00073 < 0: if INFO = -i, the i-th argument had an illegal value 00074 > 0: the algorithm failed to find all of the eigenvalues in 00075 a total of 30*N iterations; if INFO = i, then i 00076 elements of E have not converged to zero. 00077 00078 ===================================================================== 00079 00080 00081 Test the input parameters. 00082 00083 Parameter adjustments */ 00084 /* Table of constant values */ 00085 integer c__0 = 0; 00086 integer c__1 = 1; 00087 Treal c_b32 = 1.; 00088 00089 /* System generated locals */ 00090 integer i__1; 00091 Treal d__1, d__2, d__3; 00092 /* Local variables */ 00093 Treal oldc; 00094 integer lend, jtot; 00095 Treal c__; 00096 integer i__, l, m; 00097 Treal p, gamma, r__, s, alpha, sigma, anorm; 00098 integer l1; 00099 Treal bb; 00100 integer iscale; 00101 Treal oldgam, safmin; 00102 Treal safmax; 00103 integer lendsv; 00104 Treal ssfmin; 00105 integer nmaxit; 00106 Treal ssfmax, rt1, rt2, eps, rte; 00107 integer lsv; 00108 Treal eps2; 00109 00110 00111 --e; 00112 --d__; 00113 00114 /* Function Body */ 00115 *info = 0; 00116 00117 /* Quick return if possible */ 00118 00119 if (*n < 0) { 00120 *info = -1; 00121 i__1 = -(*info); 00122 template_blas_erbla("STERF ", &i__1); 00123 return 0; 00124 } 00125 if (*n <= 1) { 00126 return 0; 00127 } 00128 00129 /* Determine the unit roundoff for this environment. */ 00130 00131 eps = template_lapack_lamch("E", (Treal)0); 00132 /* Computing 2nd power */ 00133 d__1 = eps; 00134 eps2 = d__1 * d__1; 00135 safmin = template_lapack_lamch("S", (Treal)0); 00136 safmax = 1. / safmin; 00137 ssfmax = template_blas_sqrt(safmax) / 3.; 00138 ssfmin = template_blas_sqrt(safmin) / eps2; 00139 00140 /* Compute the eigenvalues of the tridiagonal matrix. */ 00141 00142 nmaxit = *n * 30; 00143 sigma = 0.; 00144 jtot = 0; 00145 00146 /* Determine where the matrix splits and choose QL or QR iteration 00147 for each block, according to whether top or bottom diagonal 00148 element is smaller. */ 00149 00150 l1 = 1; 00151 00152 L10: 00153 if (l1 > *n) { 00154 goto L170; 00155 } 00156 if (l1 > 1) { 00157 e[l1 - 1] = 0.; 00158 } 00159 i__1 = *n - 1; 00160 for (m = l1; m <= i__1; ++m) { 00161 if ((d__3 = e[m], absMACRO(d__3)) <= template_blas_sqrt((d__1 = d__[m], absMACRO(d__1))) * 00162 template_blas_sqrt((d__2 = d__[m + 1], absMACRO(d__2))) * eps) { 00163 e[m] = 0.; 00164 goto L30; 00165 } 00166 /* L20: */ 00167 } 00168 m = *n; 00169 00170 L30: 00171 l = l1; 00172 lsv = l; 00173 lend = m; 00174 lendsv = lend; 00175 l1 = m + 1; 00176 if (lend == l) { 00177 goto L10; 00178 } 00179 00180 /* Scale submatrix in rows and columns L to LEND */ 00181 00182 i__1 = lend - l + 1; 00183 anorm = template_lapack_lanst("I", &i__1, &d__[l], &e[l]); 00184 iscale = 0; 00185 if (anorm > ssfmax) { 00186 iscale = 1; 00187 i__1 = lend - l + 1; 00188 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n, 00189 info); 00190 i__1 = lend - l; 00191 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n, 00192 info); 00193 } else if (anorm < ssfmin) { 00194 iscale = 2; 00195 i__1 = lend - l + 1; 00196 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n, 00197 info); 00198 i__1 = lend - l; 00199 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n, 00200 info); 00201 } 00202 00203 i__1 = lend - 1; 00204 for (i__ = l; i__ <= i__1; ++i__) { 00205 /* Computing 2nd power */ 00206 d__1 = e[i__]; 00207 e[i__] = d__1 * d__1; 00208 /* L40: */ 00209 } 00210 00211 /* Choose between QL and QR iteration */ 00212 00213 if ((d__1 = d__[lend], absMACRO(d__1)) < (d__2 = d__[l], absMACRO(d__2))) { 00214 lend = lsv; 00215 l = lendsv; 00216 } 00217 00218 if (lend >= l) { 00219 00220 /* QL Iteration 00221 00222 Look for small subdiagonal element. */ 00223 00224 L50: 00225 if (l != lend) { 00226 i__1 = lend - 1; 00227 for (m = l; m <= i__1; ++m) { 00228 if ((d__2 = e[m], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m 00229 + 1], absMACRO(d__1))) { 00230 goto L70; 00231 } 00232 /* L60: */ 00233 } 00234 } 00235 m = lend; 00236 00237 L70: 00238 if (m < lend) { 00239 e[m] = 0.; 00240 } 00241 p = d__[l]; 00242 if (m == l) { 00243 goto L90; 00244 } 00245 00246 /* If remaining matrix is 2 by 2, use DLAE2 to compute its 00247 eigenvalues. */ 00248 00249 if (m == l + 1) { 00250 rte = template_blas_sqrt(e[l]); 00251 template_lapack_lae2(&d__[l], &rte, &d__[l + 1], &rt1, &rt2); 00252 d__[l] = rt1; 00253 d__[l + 1] = rt2; 00254 e[l] = 0.; 00255 l += 2; 00256 if (l <= lend) { 00257 goto L50; 00258 } 00259 goto L150; 00260 } 00261 00262 if (jtot == nmaxit) { 00263 goto L150; 00264 } 00265 ++jtot; 00266 00267 /* Form shift. */ 00268 00269 rte = template_blas_sqrt(e[l]); 00270 sigma = (d__[l + 1] - p) / (rte * 2.); 00271 r__ = template_lapack_lapy2(&sigma, &c_b32); 00272 sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma)); 00273 00274 c__ = 1.; 00275 s = 0.; 00276 gamma = d__[m] - sigma; 00277 p = gamma * gamma; 00278 00279 /* Inner loop */ 00280 00281 i__1 = l; 00282 for (i__ = m - 1; i__ >= i__1; --i__) { 00283 bb = e[i__]; 00284 r__ = p + bb; 00285 if (i__ != m - 1) { 00286 e[i__ + 1] = s * r__; 00287 } 00288 oldc = c__; 00289 c__ = p / r__; 00290 s = bb / r__; 00291 oldgam = gamma; 00292 alpha = d__[i__]; 00293 gamma = c__ * (alpha - sigma) - s * oldgam; 00294 d__[i__ + 1] = oldgam + (alpha - gamma); 00295 if (c__ != 0.) { 00296 p = gamma * gamma / c__; 00297 } else { 00298 p = oldc * bb; 00299 } 00300 /* L80: */ 00301 } 00302 00303 e[l] = s * p; 00304 d__[l] = sigma + gamma; 00305 goto L50; 00306 00307 /* Eigenvalue found. */ 00308 00309 L90: 00310 d__[l] = p; 00311 00312 ++l; 00313 if (l <= lend) { 00314 goto L50; 00315 } 00316 goto L150; 00317 00318 } else { 00319 00320 /* QR Iteration 00321 00322 Look for small superdiagonal element. */ 00323 00324 L100: 00325 i__1 = lend + 1; 00326 for (m = l; m >= i__1; --m) { 00327 if ((d__2 = e[m - 1], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m 00328 - 1], absMACRO(d__1))) { 00329 goto L120; 00330 } 00331 /* L110: */ 00332 } 00333 m = lend; 00334 00335 L120: 00336 if (m > lend) { 00337 e[m - 1] = 0.; 00338 } 00339 p = d__[l]; 00340 if (m == l) { 00341 goto L140; 00342 } 00343 00344 /* If remaining matrix is 2 by 2, use DLAE2 to compute its 00345 eigenvalues. */ 00346 00347 if (m == l - 1) { 00348 rte = template_blas_sqrt(e[l - 1]); 00349 template_lapack_lae2(&d__[l], &rte, &d__[l - 1], &rt1, &rt2); 00350 d__[l] = rt1; 00351 d__[l - 1] = rt2; 00352 e[l - 1] = 0.; 00353 l += -2; 00354 if (l >= lend) { 00355 goto L100; 00356 } 00357 goto L150; 00358 } 00359 00360 if (jtot == nmaxit) { 00361 goto L150; 00362 } 00363 ++jtot; 00364 00365 /* Form shift. */ 00366 00367 rte = template_blas_sqrt(e[l - 1]); 00368 sigma = (d__[l - 1] - p) / (rte * 2.); 00369 r__ = template_lapack_lapy2(&sigma, &c_b32); 00370 sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma)); 00371 00372 c__ = 1.; 00373 s = 0.; 00374 gamma = d__[m] - sigma; 00375 p = gamma * gamma; 00376 00377 /* Inner loop */ 00378 00379 i__1 = l - 1; 00380 for (i__ = m; i__ <= i__1; ++i__) { 00381 bb = e[i__]; 00382 r__ = p + bb; 00383 if (i__ != m) { 00384 e[i__ - 1] = s * r__; 00385 } 00386 oldc = c__; 00387 c__ = p / r__; 00388 s = bb / r__; 00389 oldgam = gamma; 00390 alpha = d__[i__ + 1]; 00391 gamma = c__ * (alpha - sigma) - s * oldgam; 00392 d__[i__] = oldgam + (alpha - gamma); 00393 if (c__ != 0.) { 00394 p = gamma * gamma / c__; 00395 } else { 00396 p = oldc * bb; 00397 } 00398 /* L130: */ 00399 } 00400 00401 e[l - 1] = s * p; 00402 d__[l] = sigma + gamma; 00403 goto L100; 00404 00405 /* Eigenvalue found. */ 00406 00407 L140: 00408 d__[l] = p; 00409 00410 --l; 00411 if (l >= lend) { 00412 goto L100; 00413 } 00414 goto L150; 00415 00416 } 00417 00418 /* Undo scaling if necessary */ 00419 00420 L150: 00421 if (iscale == 1) { 00422 i__1 = lendsv - lsv + 1; 00423 template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv], 00424 n, info); 00425 } 00426 if (iscale == 2) { 00427 i__1 = lendsv - lsv + 1; 00428 template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv], 00429 n, info); 00430 } 00431 00432 /* Check for no convergence to an eigenvalue after a total 00433 of N*MAXIT iterations. */ 00434 00435 if (jtot < nmaxit) { 00436 goto L10; 00437 } 00438 i__1 = *n - 1; 00439 for (i__ = 1; i__ <= i__1; ++i__) { 00440 if (e[i__] != 0.) { 00441 ++(*info); 00442 } 00443 /* L160: */ 00444 } 00445 goto L180; 00446 00447 /* Sort eigenvalues in increasing order. */ 00448 00449 L170: 00450 template_lapack_lasrt("I", n, &d__[1], info); 00451 00452 L180: 00453 return 0; 00454 00455 /* End of DSTERF */ 00456 00457 } /* dsterf_ */ 00458 00459 #endif