template_lapack_larrc.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_LARRC_HEADER
00036 #define TEMPLATE_LAPACK_LARRC_HEADER
00037 
00038 template<class Treal>
00039 int template_lapack_larrc(const char *jobt, const integer *n, const Treal *vl, 
00040         const Treal *vu, Treal *d__, Treal *e, Treal *pivmin, 
00041         integer *eigcnt, integer *lcnt, integer *rcnt, integer *info)
00042 {
00043     /* System generated locals */
00044     integer i__1;
00045     Treal d__1;
00046 
00047     /* Local variables */
00048     integer i__;
00049     Treal sl, su, tmp, tmp2;
00050     logical matt;
00051     Treal lpivot, rpivot;
00052 
00053 
00054 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00055 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00056 /*     November 2006 */
00057 
00058 /*     .. Scalar Arguments .. */
00059 /*     .. */
00060 /*     .. Array Arguments .. */
00061 /*     .. */
00062 
00063 /*  Purpose */
00064 /*  ======= */
00065 
00066 /*  Find the number of eigenvalues of the symmetric tridiagonal matrix T */
00067 /*  that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T */
00068 /*  if JOBT = 'L'. */
00069 
00070 /*  Arguments */
00071 /*  ========= */
00072 
00073 /*  JOBT    (input) CHARACTER*1 */
00074 /*          = 'T':  Compute Sturm count for matrix T. */
00075 /*          = 'L':  Compute Sturm count for matrix L D L^T. */
00076 
00077 /*  N       (input) INTEGER */
00078 /*          The order of the matrix. N > 0. */
00079 
00080 /*  VL      (input) DOUBLE PRECISION */
00081 /*  VU      (input) DOUBLE PRECISION */
00082 /*          The lower and upper bounds for the eigenvalues. */
00083 
00084 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
00085 /*          JOBT = 'T': The N diagonal elements of the tridiagonal matrix T. */
00086 /*          JOBT = 'L': The N diagonal elements of the diagonal matrix D. */
00087 
00088 /*  E       (input) DOUBLE PRECISION array, dimension (N) */
00089 /*          JOBT = 'T': The N-1 offdiagonal elements of the matrix T. */
00090 /*          JOBT = 'L': The N-1 offdiagonal elements of the matrix L. */
00091 
00092 /*  PIVMIN  (input) DOUBLE PRECISION */
00093 /*          The minimum pivot in the Sturm sequence for T. */
00094 
00095 /*  EIGCNT  (output) INTEGER */
00096 /*          The number of eigenvalues of the symmetric tridiagonal matrix T */
00097 /*          that are in the interval (VL,VU] */
00098 
00099 /*  LCNT    (output) INTEGER */
00100 /*  RCNT    (output) INTEGER */
00101 /*          The left and right negcounts of the interval. */
00102 
00103 /*  INFO    (output) INTEGER */
00104 
00105 /*  Further Details */
00106 /*  =============== */
00107 
00108 /*  Based on contributions by */
00109 /*     Beresford Parlett, University of California, Berkeley, USA */
00110 /*     Jim Demmel, University of California, Berkeley, USA */
00111 /*     Inderjit Dhillon, University of Texas, Austin, USA */
00112 /*     Osni Marques, LBNL/NERSC, USA */
00113 /*     Christof Voemel, University of California, Berkeley, USA */
00114 
00115 /*  ===================================================================== */
00116 
00117 /*     .. Parameters .. */
00118 /*     .. */
00119 /*     .. Local Scalars .. */
00120 /*     .. */
00121 /*     .. External Functions .. */
00122 /*     .. */
00123 /*     .. Executable Statements .. */
00124 
00125     /* Parameter adjustments */
00126     --e;
00127     --d__;
00128 
00129     /* Function Body */
00130     *info = 0;
00131     *lcnt = 0;
00132     *rcnt = 0;
00133     *eigcnt = 0;
00134     matt = template_blas_lsame(jobt, "T");
00135     if (matt) {
00136 /*        Sturm sequence count on T */
00137         lpivot = d__[1] - *vl;
00138         rpivot = d__[1] - *vu;
00139         if (lpivot <= 0.) {
00140             ++(*lcnt);
00141         }
00142         if (rpivot <= 0.) {
00143             ++(*rcnt);
00144         }
00145         i__1 = *n - 1;
00146         for (i__ = 1; i__ <= i__1; ++i__) {
00147 /* Computing 2nd power */
00148             d__1 = e[i__];
00149             tmp = d__1 * d__1;
00150             lpivot = d__[i__ + 1] - *vl - tmp / lpivot;
00151             rpivot = d__[i__ + 1] - *vu - tmp / rpivot;
00152             if (lpivot <= 0.) {
00153                 ++(*lcnt);
00154             }
00155             if (rpivot <= 0.) {
00156                 ++(*rcnt);
00157             }
00158 /* L10: */
00159         }
00160     } else {
00161 /*        Sturm sequence count on L D L^T */
00162         sl = -(*vl);
00163         su = -(*vu);
00164         i__1 = *n - 1;
00165         for (i__ = 1; i__ <= i__1; ++i__) {
00166             lpivot = d__[i__] + sl;
00167             rpivot = d__[i__] + su;
00168             if (lpivot <= 0.) {
00169                 ++(*lcnt);
00170             }
00171             if (rpivot <= 0.) {
00172                 ++(*rcnt);
00173             }
00174             tmp = e[i__] * d__[i__] * e[i__];
00175 
00176             tmp2 = tmp / lpivot;
00177             if (tmp2 == 0.) {
00178                 sl = tmp - *vl;
00179             } else {
00180                 sl = sl * tmp2 - *vl;
00181             }
00182 
00183             tmp2 = tmp / rpivot;
00184             if (tmp2 == 0.) {
00185                 su = tmp - *vu;
00186             } else {
00187                 su = su * tmp2 - *vu;
00188             }
00189 /* L20: */
00190         }
00191         lpivot = d__[*n] + sl;
00192         rpivot = d__[*n] + su;
00193         if (lpivot <= 0.) {
00194             ++(*lcnt);
00195         }
00196         if (rpivot <= 0.) {
00197             ++(*rcnt);
00198         }
00199     }
00200     *eigcnt = *rcnt - *lcnt;
00201     return 0;
00202 
00203 /*     end of DLARRC */
00204 
00205 } /* dlarrc_ */
00206 
00207 #endif

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1