template_lapack_syev.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_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

Generated on Mon Sep 17 14:32:56 2012 for ergo by  doxygen 1.4.7