ergo
|
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