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_SYGV_HEADER 00036 #define TEMPLATE_LAPACK_SYGV_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_sygv(const integer *itype, const char *jobz, const char *uplo, const integer * 00041 n, Treal *a, const integer *lda, Treal *b, const integer *ldb, 00042 Treal *w, Treal *work, const integer *lwork, integer *info) 00043 { 00044 00045 //printf("entering template_lapack_sygv\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 DSYGV computes all the eigenvalues, and optionally, the eigenvectors 00057 of a real generalized symmetric-definite eigenproblem, of the form 00058 A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. 00059 Here A and B are assumed to be symmetric and B is also 00060 positive definite. 00061 00062 Arguments 00063 ========= 00064 00065 ITYPE (input) INTEGER 00066 Specifies the problem type to be solved: 00067 = 1: A*x = (lambda)*B*x 00068 = 2: A*B*x = (lambda)*x 00069 = 3: B*A*x = (lambda)*x 00070 00071 JOBZ (input) CHARACTER*1 00072 = 'N': Compute eigenvalues only; 00073 = 'V': Compute eigenvalues and eigenvectors. 00074 00075 UPLO (input) CHARACTER*1 00076 = 'U': Upper triangles of A and B are stored; 00077 = 'L': Lower triangles of A and B are stored. 00078 00079 N (input) INTEGER 00080 The order of the matrices A and B. N >= 0. 00081 00082 A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00083 On entry, the symmetric matrix A. If UPLO = 'U', the 00084 leading N-by-N upper triangular part of A contains the 00085 upper triangular part of the matrix A. If UPLO = 'L', 00086 the leading N-by-N lower triangular part of A contains 00087 the lower triangular part of the matrix A. 00088 00089 On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00090 matrix Z of eigenvectors. The eigenvectors are normalized 00091 as follows: 00092 if ITYPE = 1 or 2, Z**T*B*Z = I; 00093 if ITYPE = 3, Z**T*inv(B)*Z = I. 00094 If JOBZ = 'N', then on exit the upper triangle (if UPLO='U') 00095 or the lower triangle (if UPLO='L') of A, including the 00096 diagonal, is destroyed. 00097 00098 LDA (input) INTEGER 00099 The leading dimension of the array A. LDA >= max(1,N). 00100 00101 B (input/output) DOUBLE PRECISION array, dimension (LDB, N) 00102 On entry, the symmetric positive definite matrix B. 00103 If UPLO = 'U', the leading N-by-N upper triangular part of B 00104 contains the upper triangular part of the matrix B. 00105 If UPLO = 'L', the leading N-by-N lower triangular part of B 00106 contains the lower triangular part of the matrix B. 00107 00108 On exit, if INFO <= N, the part of B containing the matrix is 00109 overwritten by the triangular factor U or L from the Cholesky 00110 factorization B = U**T*U or B = L*L**T. 00111 00112 LDB (input) INTEGER 00113 The leading dimension of the array B. LDB >= max(1,N). 00114 00115 W (output) DOUBLE PRECISION array, dimension (N) 00116 If INFO = 0, the eigenvalues in ascending order. 00117 00118 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 00119 On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00120 00121 LWORK (input) INTEGER 00122 The length of the array WORK. LWORK >= max(1,3*N-1). 00123 For optimal efficiency, LWORK >= (NB+2)*N, 00124 where NB is the blocksize for DSYTRD returned by ILAENV. 00125 00126 If LWORK = -1, then a workspace query is assumed; the routine 00127 only calculates the optimal size of the WORK array, returns 00128 this value as the first entry of the WORK array, and no error 00129 message related to LWORK is issued by XERBLA. 00130 00131 INFO (output) INTEGER 00132 = 0: successful exit 00133 < 0: if INFO = -i, the i-th argument had an illegal value 00134 > 0: DPOTRF or DSYEV returned an error code: 00135 <= N: if INFO = i, DSYEV failed to converge; 00136 i off-diagonal elements of an intermediate 00137 tridiagonal form did not converge to zero; 00138 > N: if INFO = N + i, for 1 <= i <= N, then the leading 00139 minor of order i of B is not positive definite. 00140 The factorization of B could not be completed and 00141 no eigenvalues or eigenvectors were computed. 00142 00143 ===================================================================== 00144 00145 00146 Test the input parameters. 00147 00148 Parameter adjustments */ 00149 /* Table of constant values */ 00150 integer c__1 = 1; 00151 integer c_n1 = -1; 00152 Treal c_b16 = 1.; 00153 00154 /* System generated locals */ 00155 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2; 00156 /* Local variables */ 00157 integer neig; 00158 char trans[1]; 00159 logical upper; 00160 logical wantz; 00161 integer nb; 00162 integer lwkopt; 00163 logical lquery; 00164 00165 00166 a_dim1 = *lda; 00167 a_offset = 1 + a_dim1 * 1; 00168 a -= a_offset; 00169 b_dim1 = *ldb; 00170 b_offset = 1 + b_dim1 * 1; 00171 b -= b_offset; 00172 --w; 00173 --work; 00174 00175 /* Initialization added by Elias to get rid of compiler warnings. */ 00176 lwkopt = 0; 00177 /* Function Body */ 00178 wantz = template_blas_lsame(jobz, "V"); 00179 upper = template_blas_lsame(uplo, "U"); 00180 lquery = *lwork == -1; 00181 00182 *info = 0; 00183 if (*itype < 1 || *itype > 3) { 00184 *info = -1; 00185 } else if (! (wantz || template_blas_lsame(jobz, "N"))) { 00186 *info = -2; 00187 } else if (! (upper || template_blas_lsame(uplo, "L"))) { 00188 *info = -3; 00189 } else if (*n < 0) { 00190 *info = -4; 00191 } else if (*lda < maxMACRO(1,*n)) { 00192 *info = -6; 00193 } else if (*ldb < maxMACRO(1,*n)) { 00194 *info = -8; 00195 } else /* if(complicated condition) */ { 00196 /* Computing MAX */ 00197 i__1 = 1, i__2 = *n * 3 - 1; 00198 if (*lwork < maxMACRO(i__1,i__2) && ! lquery) { 00199 *info = -11; 00200 } 00201 } 00202 00203 if (*info == 0) { 00204 nb = template_lapack_ilaenv(&c__1, "DSYTRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, 00205 (ftnlen)1); 00206 lwkopt = (nb + 2) * *n; 00207 work[1] = (Treal) lwkopt; 00208 } 00209 00210 if (*info != 0) { 00211 i__1 = -(*info); 00212 template_blas_erbla("SYGV ", &i__1); 00213 return 0; 00214 } else if (lquery) { 00215 return 0; 00216 } 00217 00218 /* Quick return if possible */ 00219 00220 if (*n == 0) { 00221 return 0; 00222 } 00223 00224 /* Form a Cholesky factorization of B. */ 00225 00226 //printf("calling template_lapack_potrf\n"); 00227 template_lapack_potrf(uplo, n, &b[b_offset], ldb, info); 00228 //printf("template_lapack_potrf returned\n"); 00229 00230 00231 if (*info != 0) { 00232 *info = *n + *info; 00233 return 0; 00234 } 00235 00236 /* Transform problem to standard eigenvalue problem and solve. */ 00237 00238 //printf("calling template_lapack_sygst\n"); 00239 template_lapack_sygst(itype, uplo, n, &a[a_offset], lda, &b[b_offset], ldb, info); 00240 //printf("template_lapack_sygst returned\n"); 00241 00242 //printf("calling template_lapack_syev\n"); 00243 template_lapack_syev(jobz, uplo, n, &a[a_offset], lda, &w[1], &work[1], lwork, info); 00244 //printf("template_lapack_syev returned\n"); 00245 00246 if (wantz) { 00247 00248 /* Backtransform eigenvectors to the original problem. */ 00249 00250 neig = *n; 00251 if (*info > 0) { 00252 neig = *info - 1; 00253 } 00254 if (*itype == 1 || *itype == 2) { 00255 00256 /* For A*x=(lambda)*B*x and A*B*x=(lambda)*x; 00257 backtransform eigenvectors: x = inv(L)'*y or inv(U)*y */ 00258 00259 if (upper) { 00260 *(unsigned char *)trans = 'N'; 00261 } else { 00262 *(unsigned char *)trans = 'T'; 00263 } 00264 00265 template_blas_trsm("Left", uplo, trans, "Non-unit", n, &neig, &c_b16, &b[ 00266 b_offset], ldb, &a[a_offset], lda); 00267 00268 } else if (*itype == 3) { 00269 00270 /* For B*A*x=(lambda)*x; 00271 backtransform eigenvectors: x = L*y or U'*y */ 00272 00273 if (upper) { 00274 *(unsigned char *)trans = 'T'; 00275 } else { 00276 *(unsigned char *)trans = 'N'; 00277 } 00278 00279 template_blas_trmm("Left", uplo, trans, "Non-unit", n, &neig, &c_b16, &b[ 00280 b_offset], ldb, &a[a_offset], lda); 00281 } 00282 } 00283 00284 work[1] = (Treal) lwkopt; 00285 return 0; 00286 00287 /* End of DSYGV */ 00288 00289 } /* dsygv_ */ 00290 00291 #endif