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_LARTG_HEADER 00036 #define TEMPLATE_LAPACK_LARTG_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, 00041 Treal *sn, Treal *r__) 00042 { 00043 /* -- LAPACK auxiliary routine (version 3.0) -- 00044 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00045 Courant Institute, Argonne National Lab, and Rice University 00046 September 30, 1994 00047 00048 00049 Purpose 00050 ======= 00051 00052 DLARTG generate a plane rotation so that 00053 00054 [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. 00055 [ -SN CS ] [ G ] [ 0 ] 00056 00057 This is a slower, more accurate version of the BLAS1 routine DROTG, 00058 with the following other differences: 00059 F and G are unchanged on return. 00060 If G=0, then CS=1 and SN=0. 00061 If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any 00062 floating point operations (saves work in DBDSQR when 00063 there are zeros on the diagonal). 00064 00065 If F exceeds G in magnitude, CS will be positive. 00066 00067 Arguments 00068 ========= 00069 00070 F (input) DOUBLE PRECISION 00071 The first component of vector to be rotated. 00072 00073 G (input) DOUBLE PRECISION 00074 The second component of vector to be rotated. 00075 00076 CS (output) DOUBLE PRECISION 00077 The cosine of the rotation. 00078 00079 SN (output) DOUBLE PRECISION 00080 The sine of the rotation. 00081 00082 R (output) DOUBLE PRECISION 00083 The nonzero component of the rotated vector. 00084 00085 ===================================================================== */ 00086 /* Initialized data */ 00087 logical first = TRUE_; 00088 /* System generated locals */ 00089 integer i__1; 00090 Treal d__1, d__2; 00091 /* Local variables */ 00092 integer i__; 00093 Treal scale; 00094 integer count; 00095 Treal f1, g1, safmn2, safmx2; 00096 Treal safmin, eps; 00097 00098 00099 00100 if (first) { 00101 first = FALSE_; 00102 safmin = template_lapack_lamch("S", (Treal)0); 00103 eps = template_lapack_lamch("E", (Treal)0); 00104 d__1 = template_lapack_lamch("B", (Treal)0); 00105 i__1 = (integer) (template_blas_log(safmin / eps) / template_blas_log(template_lapack_lamch("B", (Treal)0)) / 00106 2.); 00107 safmn2 = template_lapack_pow_di(&d__1, &i__1); 00108 safmx2 = 1. / safmn2; 00109 } 00110 if (*g == 0.) { 00111 *cs = 1.; 00112 *sn = 0.; 00113 *r__ = *f; 00114 } else if (*f == 0.) { 00115 *cs = 0.; 00116 *sn = 1.; 00117 *r__ = *g; 00118 } else { 00119 f1 = *f; 00120 g1 = *g; 00121 /* Computing MAX */ 00122 d__1 = absMACRO(f1), d__2 = absMACRO(g1); 00123 scale = maxMACRO(d__1,d__2); 00124 if (scale >= safmx2) { 00125 count = 0; 00126 L10: 00127 ++count; 00128 f1 *= safmn2; 00129 g1 *= safmn2; 00130 /* Computing MAX */ 00131 d__1 = absMACRO(f1), d__2 = absMACRO(g1); 00132 scale = maxMACRO(d__1,d__2); 00133 if (scale >= safmx2) { 00134 goto L10; 00135 } 00136 /* Computing 2nd power */ 00137 d__1 = f1; 00138 /* Computing 2nd power */ 00139 d__2 = g1; 00140 *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2); 00141 *cs = f1 / *r__; 00142 *sn = g1 / *r__; 00143 i__1 = count; 00144 for (i__ = 1; i__ <= i__1; ++i__) { 00145 *r__ *= safmx2; 00146 /* L20: */ 00147 } 00148 } else if (scale <= safmn2) { 00149 count = 0; 00150 L30: 00151 ++count; 00152 f1 *= safmx2; 00153 g1 *= safmx2; 00154 /* Computing MAX */ 00155 d__1 = absMACRO(f1), d__2 = absMACRO(g1); 00156 scale = maxMACRO(d__1,d__2); 00157 if (scale <= safmn2) { 00158 goto L30; 00159 } 00160 /* Computing 2nd power */ 00161 d__1 = f1; 00162 /* Computing 2nd power */ 00163 d__2 = g1; 00164 *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2); 00165 *cs = f1 / *r__; 00166 *sn = g1 / *r__; 00167 i__1 = count; 00168 for (i__ = 1; i__ <= i__1; ++i__) { 00169 *r__ *= safmn2; 00170 /* L40: */ 00171 } 00172 } else { 00173 /* Computing 2nd power */ 00174 d__1 = f1; 00175 /* Computing 2nd power */ 00176 d__2 = g1; 00177 *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2); 00178 *cs = f1 / *r__; 00179 *sn = g1 / *r__; 00180 } 00181 if (absMACRO(*f) > absMACRO(*g) && *cs < 0.) { 00182 *cs = -(*cs); 00183 *sn = -(*sn); 00184 *r__ = -(*r__); 00185 } 00186 } 00187 return 0; 00188 00189 /* End of DLARTG */ 00190 00191 } /* dlartg_ */ 00192 00193 #endif