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_SYEV_HEADER 00036 #define TEMPLATE_LAPACK_SYEV_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_syev(const char *jobz, const char *uplo, const integer *n, Treal *a, 00041 const integer *lda, Treal *w, Treal *work, const integer *lwork, 00042 integer *info) 00043 { 00044 00045 //printf("entering template_lapack_syev\n"); 00046 00047 /* -- LAPACK driver routine (version 3.0) -- 00048 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00049 Courant Institute, Argonne National Lab, and Rice University 00050 June 30, 1999 00051 00052 00053 Purpose 00054 ======= 00055 00056 DSYEV computes all eigenvalues and, optionally, eigenvectors of a 00057 real symmetric matrix A. 00058 00059 Arguments 00060 ========= 00061 00062 JOBZ (input) CHARACTER*1 00063 = 'N': Compute eigenvalues only; 00064 = 'V': Compute eigenvalues and eigenvectors. 00065 00066 UPLO (input) CHARACTER*1 00067 = 'U': Upper triangle of A is stored; 00068 = 'L': Lower triangle of A is stored. 00069 00070 N (input) INTEGER 00071 The order of the matrix A. N >= 0. 00072 00073 A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00074 On entry, the symmetric matrix A. If UPLO = 'U', the 00075 leading N-by-N upper triangular part of A contains the 00076 upper triangular part of the matrix A. If UPLO = 'L', 00077 the leading N-by-N lower triangular part of A contains 00078 the lower triangular part of the matrix A. 00079 On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00080 orthonormal eigenvectors of the matrix A. 00081 If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 00082 or the upper triangle (if UPLO='U') of A, including the 00083 diagonal, is destroyed. 00084 00085 LDA (input) INTEGER 00086 The leading dimension of the array A. LDA >= max(1,N). 00087 00088 W (output) DOUBLE PRECISION array, dimension (N) 00089 If INFO = 0, the eigenvalues in ascending order. 00090 00091 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 00092 On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00093 00094 LWORK (input) INTEGER 00095 The length of the array WORK. LWORK >= max(1,3*N-1). 00096 For optimal efficiency, LWORK >= (NB+2)*N, 00097 where NB is the blocksize for DSYTRD returned by ILAENV. 00098 00099 If LWORK = -1, then a workspace query is assumed; the routine 00100 only calculates the optimal size of the WORK array, returns 00101 this value as the first entry of the WORK array, and no error 00102 message related to LWORK is issued by XERBLA. 00103 00104 INFO (output) INTEGER 00105 = 0: successful exit 00106 < 0: if INFO = -i, the i-th argument had an illegal value 00107 > 0: if INFO = i, the algorithm failed to converge; i 00108 off-diagonal elements of an intermediate tridiagonal 00109 form did not converge to zero. 00110 00111 ===================================================================== 00112 00113 00114 Test the input parameters. 00115 00116 Parameter adjustments */ 00117 /* Table of constant values */ 00118 integer c__1 = 1; 00119 integer c_n1 = -1; 00120 integer c__0 = 0; 00121 Treal c_b17 = 1.; 00122 00123 /* System generated locals */ 00124 integer a_dim1, a_offset, i__1, i__2; 00125 Treal d__1; 00126 /* Local variables */ 00127 integer inde; 00128 Treal anrm; 00129 integer imax; 00130 Treal rmin, rmax; 00131 Treal sigma; 00132 integer iinfo; 00133 logical lower, wantz; 00134 integer nb; 00135 integer iscale; 00136 Treal safmin; 00137 Treal bignum; 00138 integer indtau; 00139 integer indwrk; 00140 integer llwork; 00141 Treal smlnum; 00142 integer lwkopt; 00143 logical lquery; 00144 Treal eps; 00145 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00146 00147 00148 a_dim1 = *lda; 00149 a_offset = 1 + a_dim1 * 1; 00150 a -= a_offset; 00151 --w; 00152 --work; 00153 00154 /* Initialization added by Elias to get rid of compiler warnings. */ 00155 lwkopt = 0; 00156 /* Function Body */ 00157 wantz = template_blas_lsame(jobz, "V"); 00158 lower = template_blas_lsame(uplo, "L"); 00159 lquery = *lwork == -1; 00160 00161 *info = 0; 00162 if (! (wantz || template_blas_lsame(jobz, "N"))) { 00163 *info = -1; 00164 } else if (! (lower || template_blas_lsame(uplo, "U"))) { 00165 *info = -2; 00166 } else if (*n < 0) { 00167 *info = -3; 00168 } else if (*lda < maxMACRO(1,*n)) { 00169 *info = -5; 00170 } else /* if(complicated condition) */ { 00171 /* Computing MAX */ 00172 i__1 = 1, i__2 = *n * 3 - 1; 00173 if (*lwork < maxMACRO(i__1,i__2) && ! lquery) { 00174 *info = -8; 00175 } 00176 } 00177 00178 if (*info == 0) { 00179 nb = template_lapack_ilaenv(&c__1, "DSYTRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, 00180 (ftnlen)1); 00181 /* Computing MAX */ 00182 i__1 = 1, i__2 = (nb + 2) * *n; 00183 lwkopt = maxMACRO(i__1,i__2); 00184 work[1] = (Treal) lwkopt; 00185 } 00186 00187 if (*info != 0) { 00188 i__1 = -(*info); 00189 template_blas_erbla("SYEV ", &i__1); 00190 return 0; 00191 } else if (lquery) { 00192 return 0; 00193 } 00194 00195 /* Quick return if possible */ 00196 00197 if (*n == 0) { 00198 work[1] = 1.; 00199 return 0; 00200 } 00201 00202 if (*n == 1) { 00203 w[1] = a_ref(1, 1); 00204 work[1] = 3.; 00205 if (wantz) { 00206 a_ref(1, 1) = 1.; 00207 } 00208 return 0; 00209 } 00210 00211 /* Get machine constants. */ 00212 00213 //printf("before getting machine constants.\n"); 00214 00215 //printf("calling template_lapack_lamch for Safe minimum\n"); 00216 safmin = template_lapack_lamch("Safe minimum", (Treal)0); 00217 //printf("template_lapack_lamch for Safe minimum returned\n"); 00218 00219 eps = template_lapack_lamch("Precision", (Treal)0); 00220 00221 //printf("after getting machine constants.\n"); 00222 00223 smlnum = safmin / eps; 00224 bignum = 1. / smlnum; 00225 rmin = template_blas_sqrt(smlnum); 00226 rmax = template_blas_sqrt(bignum); 00227 00228 /* Scale matrix to allowable range, if necessary. */ 00229 00230 anrm = template_lapack_lansy("M", uplo, n, &a[a_offset], lda, &work[1]); 00231 iscale = 0; 00232 if (anrm > 0. && anrm < rmin) { 00233 iscale = 1; 00234 sigma = rmin / anrm; 00235 } else if (anrm > rmax) { 00236 iscale = 1; 00237 sigma = rmax / anrm; 00238 } 00239 if (iscale == 1) { 00240 template_lapack_lascl(uplo, &c__0, &c__0, &c_b17, &sigma, n, n, &a[a_offset], lda, 00241 info); 00242 } 00243 00244 /* Call DSYTRD to reduce symmetric matrix to tridiagonal form. */ 00245 00246 inde = 1; 00247 indtau = inde + *n; 00248 indwrk = indtau + *n; 00249 llwork = *lwork - indwrk + 1; 00250 template_lapack_sytrd(uplo, n, &a[a_offset], lda, &w[1], &work[inde], &work[indtau], & 00251 work[indwrk], &llwork, &iinfo); 00252 00253 /* For eigenvalues only, call DSTERF. For eigenvectors, first call 00254 DORGTR to generate the orthogonal matrix, then call DSTEQR. */ 00255 00256 if (! wantz) { 00257 template_lapack_sterf(n, &w[1], &work[inde], info); 00258 } else { 00259 template_lapack_orgtr(uplo, n, &a[a_offset], lda, &work[indtau], &work[indwrk], & 00260 llwork, &iinfo); 00261 template_lapack_steqr(jobz, n, &w[1], &work[inde], &a[a_offset], lda, &work[indtau], 00262 info); 00263 } 00264 00265 /* If matrix was scaled, then rescale eigenvalues appropriately. */ 00266 00267 if (iscale == 1) { 00268 if (*info == 0) { 00269 imax = *n; 00270 } else { 00271 imax = *info - 1; 00272 } 00273 d__1 = 1. / sigma; 00274 template_blas_scal(&imax, &d__1, &w[1], &c__1); 00275 } 00276 00277 /* Set WORK(1) to optimal workspace size. */ 00278 00279 work[1] = (Treal) lwkopt; 00280 00281 return 0; 00282 00283 /* End of DSYEV */ 00284 00285 } /* dsyev_ */ 00286 00287 #undef a_ref 00288 00289 00290 #endif