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