template_lapack_lasq6.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_LASQ6_HEADER
00036 #define TEMPLATE_LAPACK_LASQ6_HEADER
00037 
00038 template<class Treal>
00039 int template_lapack_lasq6(integer *i0, integer *n0, Treal *z__, 
00040         integer *pp, Treal *dmin__, Treal *dmin1, Treal *dmin2, 
00041          Treal *dn, Treal *dnm1, Treal *dnm2)
00042 {
00043     /* System generated locals */
00044     integer i__1;
00045     Treal d__1, d__2;
00046 
00047     /* Local variables */
00048     Treal d__;
00049     integer j4, j4p2;
00050     Treal emin, temp;
00051     Treal safmin;
00052 
00053 
00054 /*  -- LAPACK routine (version 3.2)                                    -- */
00055 
00056 /*  -- Contributed by Osni Marques of the Lawrence Berkeley National   -- */
00057 /*  -- Laboratory and Beresford Parlett of the Univ. of California at  -- */
00058 /*  -- Berkeley                                                        -- */
00059 /*  -- November 2008                                                   -- */
00060 
00061 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00062 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00063 
00064 /*     .. Scalar Arguments .. */
00065 /*     .. */
00066 /*     .. Array Arguments .. */
00067 /*     .. */
00068 
00069 /*  Purpose */
00070 /*  ======= */
00071 
00072 /*  DLASQ6 computes one dqd (shift equal to zero) transform in */
00073 /*  ping-pong form, with protection against underflow and overflow. */
00074 
00075 /*  Arguments */
00076 /*  ========= */
00077 
00078 /*  I0    (input) INTEGER */
00079 /*        First index. */
00080 
00081 /*  N0    (input) INTEGER */
00082 /*        Last index. */
00083 
00084 /*  Z     (input) DOUBLE PRECISION array, dimension ( 4*N ) */
00085 /*        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid */
00086 /*        an extra argument. */
00087 
00088 /*  PP    (input) INTEGER */
00089 /*        PP=0 for ping, PP=1 for pong. */
00090 
00091 /*  DMIN  (output) DOUBLE PRECISION */
00092 /*        Minimum value of d. */
00093 
00094 /*  DMIN1 (output) DOUBLE PRECISION */
00095 /*        Minimum value of d, excluding D( N0 ). */
00096 
00097 /*  DMIN2 (output) DOUBLE PRECISION */
00098 /*        Minimum value of d, excluding D( N0 ) and D( N0-1 ). */
00099 
00100 /*  DN    (output) DOUBLE PRECISION */
00101 /*        d(N0), the last value of d. */
00102 
00103 /*  DNM1  (output) DOUBLE PRECISION */
00104 /*        d(N0-1). */
00105 
00106 /*  DNM2  (output) DOUBLE PRECISION */
00107 /*        d(N0-2). */
00108 
00109 /*  ===================================================================== */
00110 
00111 /*     .. Parameter .. */
00112 /*     .. */
00113 /*     .. Local Scalars .. */
00114 /*     .. */
00115 /*     .. External Function .. */
00116 /*     .. */
00117 /*     .. Intrinsic Functions .. */
00118 /*     .. */
00119 /*     .. Executable Statements .. */
00120 
00121     /* Parameter adjustments */
00122     --z__;
00123 
00124     /* Function Body */
00125     if (*n0 - *i0 - 1 <= 0) {
00126         return 0;
00127     }
00128 
00129     safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00130     j4 = (*i0 << 2) + *pp - 3;
00131     emin = z__[j4 + 4];
00132     d__ = z__[j4];
00133     *dmin__ = d__;
00134 
00135     if (*pp == 0) {
00136       i__1 = (*n0 - 3) << 2;
00137         for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
00138             z__[j4 - 2] = d__ + z__[j4 - 1];
00139             if (z__[j4 - 2] == 0.) {
00140                 z__[j4] = 0.;
00141                 d__ = z__[j4 + 1];
00142                 *dmin__ = d__;
00143                 emin = 0.;
00144             } else if (safmin * z__[j4 + 1] < z__[j4 - 2] && safmin * z__[j4 
00145                     - 2] < z__[j4 + 1]) {
00146                 temp = z__[j4 + 1] / z__[j4 - 2];
00147                 z__[j4] = z__[j4 - 1] * temp;
00148                 d__ *= temp;
00149             } else {
00150                 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
00151                 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]);
00152             }
00153             *dmin__ = minMACRO(*dmin__,d__);
00154 /* Computing MIN */
00155             d__1 = emin, d__2 = z__[j4];
00156             emin = minMACRO(d__1,d__2);
00157 /* L10: */
00158         }
00159     } else {
00160       i__1 = ( *n0 - 3 ) << 2;
00161         for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
00162             z__[j4 - 3] = d__ + z__[j4];
00163             if (z__[j4 - 3] == 0.) {
00164                 z__[j4 - 1] = 0.;
00165                 d__ = z__[j4 + 2];
00166                 *dmin__ = d__;
00167                 emin = 0.;
00168             } else if (safmin * z__[j4 + 2] < z__[j4 - 3] && safmin * z__[j4 
00169                     - 3] < z__[j4 + 2]) {
00170                 temp = z__[j4 + 2] / z__[j4 - 3];
00171                 z__[j4 - 1] = z__[j4] * temp;
00172                 d__ *= temp;
00173             } else {
00174                 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
00175                 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]);
00176             }
00177             *dmin__ = minMACRO(*dmin__,d__);
00178 /* Computing MIN */
00179             d__1 = emin, d__2 = z__[j4 - 1];
00180             emin = minMACRO(d__1,d__2);
00181 /* L20: */
00182         }
00183     }
00184 
00185 /*     Unroll last two steps. */
00186 
00187     *dnm2 = d__;
00188     *dmin2 = *dmin__;
00189     j4 = ( ( *n0 - 2 ) << 2) - *pp;
00190     j4p2 = j4 + (*pp << 1) - 1;
00191     z__[j4 - 2] = *dnm2 + z__[j4p2];
00192     if (z__[j4 - 2] == 0.) {
00193         z__[j4] = 0.;
00194         *dnm1 = z__[j4p2 + 2];
00195         *dmin__ = *dnm1;
00196         emin = 0.;
00197     } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] < 
00198             z__[j4p2 + 2]) {
00199         temp = z__[j4p2 + 2] / z__[j4 - 2];
00200         z__[j4] = z__[j4p2] * temp;
00201         *dnm1 = *dnm2 * temp;
00202     } else {
00203         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
00204         *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]);
00205     }
00206     *dmin__ = minMACRO(*dmin__,*dnm1);
00207 
00208     *dmin1 = *dmin__;
00209     j4 += 4;
00210     j4p2 = j4 + (*pp << 1) - 1;
00211     z__[j4 - 2] = *dnm1 + z__[j4p2];
00212     if (z__[j4 - 2] == 0.) {
00213         z__[j4] = 0.;
00214         *dn = z__[j4p2 + 2];
00215         *dmin__ = *dn;
00216         emin = 0.;
00217     } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] < 
00218             z__[j4p2 + 2]) {
00219         temp = z__[j4p2 + 2] / z__[j4 - 2];
00220         z__[j4] = z__[j4p2] * temp;
00221         *dn = *dnm1 * temp;
00222     } else {
00223         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
00224         *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]);
00225     }
00226     *dmin__ = minMACRO(*dmin__,*dn);
00227 
00228     z__[j4 + 2] = *dn;
00229     z__[(*n0 << 2) - *pp] = emin;
00230     return 0;
00231 
00232 /*     End of DLASQ6 */
00233 
00234 } /* dlasq6_ */
00235 
00236 #endif

Generated on Mon Sep 17 14:30:40 2012 for ergo by  doxygen 1.4.7