template_lapack_larrr.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_LARRR_HEADER
00036 #define TEMPLATE_LAPACK_LARRR_HEADER
00037 
00038 template<class Treal>
00039 int template_lapack_larrr(const integer *n, Treal *d__, Treal *e, 
00040         integer *info)
00041 {
00042     /* System generated locals */
00043     integer i__1;
00044     Treal d__1;
00045 
00046 
00047     /* Local variables */
00048     integer i__;
00049     Treal eps, tmp, tmp2, rmin;
00050     Treal offdig, safmin;
00051     logical yesrel;
00052     Treal smlnum, offdig2;
00053 
00054 
00055 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00056 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00057 /*     November 2006 */
00058 
00059 /*     .. Scalar Arguments .. */
00060 /*     .. */
00061 /*     .. Array Arguments .. */
00062 /*     .. */
00063 
00064 
00065 /*  Purpose */
00066 /*  ======= */
00067 
00068 /*  Perform tests to decide whether the symmetric tridiagonal matrix T */
00069 /*  warrants expensive computations which guarantee high relative accuracy */
00070 /*  in the eigenvalues. */
00071 
00072 /*  Arguments */
00073 /*  ========= */
00074 
00075 /*  N       (input) INTEGER */
00076 /*          The order of the matrix. N > 0. */
00077 
00078 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
00079 /*          The N diagonal elements of the tridiagonal matrix T. */
00080 
00081 /*  E       (input/output) DOUBLE PRECISION array, dimension (N) */
00082 /*          On entry, the first (N-1) entries contain the subdiagonal */
00083 /*          elements of the tridiagonal matrix T; E(N) is set to ZERO. */
00084 
00085 /*  INFO    (output) INTEGER */
00086 /*          INFO = 0(default) : the matrix warrants computations preserving */
00087 /*                              relative accuracy. */
00088 /*          INFO = 1          : the matrix warrants computations guaranteeing */
00089 /*                              only absolute accuracy. */
00090 
00091 /*  Further Details */
00092 /*  =============== */
00093 
00094 /*  Based on contributions by */
00095 /*     Beresford Parlett, University of California, Berkeley, USA */
00096 /*     Jim Demmel, University of California, Berkeley, USA */
00097 /*     Inderjit Dhillon, University of Texas, Austin, USA */
00098 /*     Osni Marques, LBNL/NERSC, USA */
00099 /*     Christof Voemel, University of California, Berkeley, USA */
00100 
00101 /*  ===================================================================== */
00102 
00103 /*     .. Parameters .. */
00104 /*     .. */
00105 /*     .. Local Scalars .. */
00106 /*     .. */
00107 /*     .. External Functions .. */
00108 /*     .. */
00109 /*     .. Intrinsic Functions .. */
00110 /*     .. */
00111 /*     .. Executable Statements .. */
00112 
00113 /*     As a default, do NOT go for relative-accuracy preserving computations. */
00114     /* Parameter adjustments */
00115     --e;
00116     --d__;
00117 
00118     /* Function Body */
00119     *info = 1;
00120     safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00121     eps = template_lapack_lamch("Precision", (Treal)0);
00122     smlnum = safmin / eps;
00123     rmin = template_blas_sqrt(smlnum);
00124 /*     Tests for relative accuracy */
00125 
00126 /*     Test for scaled diagonal dominance */
00127 /*     Scale the diagonal entries to one and check whether the sum of the */
00128 /*     off-diagonals is less than one */
00129 
00130 /*     The sdd relative error bounds have a 1/(1- 2*x) factor in them, */
00131 /*     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative */
00132 /*     accuracy is promised.  In the notation of the code fragment below, */
00133 /*     1/(1 - (OFFDIG + OFFDIG2)) is the condition number. */
00134 /*     We don't think it is worth going into "sdd mode" unless the relative */
00135 /*     condition number is reasonable, not 1/macheps. */
00136 /*     The threshold should be compatible with other thresholds used in the */
00137 /*     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds */
00138 /*     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000 */
00139 /*     instead of the current OFFDIG + OFFDIG2 < 1 */
00140 
00141     yesrel = TRUE_;
00142     offdig = 0.;
00143     tmp = template_blas_sqrt((absMACRO(d__[1])));
00144     if (tmp < rmin) {
00145         yesrel = FALSE_;
00146     }
00147     if (! yesrel) {
00148         goto L11;
00149     }
00150     i__1 = *n;
00151     for (i__ = 2; i__ <= i__1; ++i__) {
00152         tmp2 = template_blas_sqrt((d__1 = d__[i__], absMACRO(d__1)));
00153         if (tmp2 < rmin) {
00154             yesrel = FALSE_;
00155         }
00156         if (! yesrel) {
00157             goto L11;
00158         }
00159         offdig2 = (d__1 = e[i__ - 1], absMACRO(d__1)) / (tmp * tmp2);
00160         if (offdig + offdig2 >= .999) {
00161             yesrel = FALSE_;
00162         }
00163         if (! yesrel) {
00164             goto L11;
00165         }
00166         tmp = tmp2;
00167         offdig = offdig2;
00168 /* L10: */
00169     }
00170 L11:
00171     if (yesrel) {
00172         *info = 0;
00173         return 0;
00174     } else {
00175     }
00176 
00177 
00178 /*     *** MORE TO BE IMPLEMENTED *** */
00179 
00180 
00181 /*     Test if the lower bidiagonal matrix L from T = L D L^T */
00182 /*     (zero shift facto) is well conditioned */
00183 
00184 
00185 /*     Test if the upper bidiagonal matrix U from T = U D U^T */
00186 /*     (zero shift facto) is well conditioned. */
00187 /*     In this case, the matrix needs to be flipped and, at the end */
00188 /*     of the eigenvector computation, the flip needs to be applied */
00189 /*     to the computed eigenvectors (and the support) */
00190 
00191 
00192     return 0;
00193 
00194 /*     END OF DLARRR */
00195 
00196 } /* dlarrr_ */
00197 
00198 #endif

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