template_lapack_sygv.h

Go to the documentation of this file.
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

Generated on Wed Nov 21 09:32:01 2012 for ergo by  doxygen 1.4.7