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_LARFG_HEADER 00036 #define TEMPLATE_LAPACK_LARFG_HEADER 00037 00038 #include "template_lapack_lapy2.h" 00039 00040 template<class Treal> 00041 int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, 00042 const integer *incx, Treal *tau) 00043 { 00044 /* -- LAPACK auxiliary routine (version 3.0) -- 00045 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00046 Courant Institute, Argonne National Lab, and Rice University 00047 September 30, 1994 00048 00049 00050 Purpose 00051 ======= 00052 00053 DLARFG generates a real elementary reflector H of order n, such 00054 that 00055 00056 H * ( alpha ) = ( beta ), H' * H = I. 00057 ( x ) ( 0 ) 00058 00059 where alpha and beta are scalars, and x is an (n-1)-element real 00060 vector. H is represented in the form 00061 00062 H = I - tau * ( 1 ) * ( 1 v' ) , 00063 ( v ) 00064 00065 where tau is a real scalar and v is a real (n-1)-element 00066 vector. 00067 00068 If the elements of x are all zero, then tau = 0 and H is taken to be 00069 the unit matrix. 00070 00071 Otherwise 1 <= tau <= 2. 00072 00073 Arguments 00074 ========= 00075 00076 N (input) INTEGER 00077 The order of the elementary reflector. 00078 00079 ALPHA (input/output) DOUBLE PRECISION 00080 On entry, the value alpha. 00081 On exit, it is overwritten with the value beta. 00082 00083 X (input/output) DOUBLE PRECISION array, dimension 00084 (1+(N-2)*abs(INCX)) 00085 On entry, the vector x. 00086 On exit, it is overwritten with the vector v. 00087 00088 INCX (input) INTEGER 00089 The increment between elements of X. INCX > 0. 00090 00091 TAU (output) DOUBLE PRECISION 00092 The value tau. 00093 00094 ===================================================================== 00095 00096 00097 Parameter adjustments */ 00098 /* System generated locals */ 00099 integer i__1; 00100 Treal d__1; 00101 /* Local variables */ 00102 Treal beta; 00103 integer j; 00104 Treal xnorm; 00105 Treal safmin, rsafmn; 00106 integer knt; 00107 00108 --x; 00109 00110 /* Function Body */ 00111 if (*n <= 1) { 00112 *tau = 0.; 00113 return 0; 00114 } 00115 00116 i__1 = *n - 1; 00117 xnorm = template_blas_nrm2(&i__1, &x[1], incx); 00118 00119 if (xnorm == 0.) { 00120 00121 /* H = I */ 00122 00123 *tau = 0.; 00124 } else { 00125 00126 /* general case */ 00127 00128 d__1 = template_lapack_lapy2(alpha, &xnorm); 00129 beta = -template_lapack_d_sign(&d__1, alpha); 00130 safmin = template_lapack_lamch("S", (Treal)0) / template_lapack_lamch("E", (Treal)0); 00131 if (absMACRO(beta) < safmin) { 00132 00133 /* XNORM, BETA may be inaccurate; scale X and recompute them */ 00134 00135 rsafmn = 1. / safmin; 00136 knt = 0; 00137 L10: 00138 ++knt; 00139 i__1 = *n - 1; 00140 template_blas_scal(&i__1, &rsafmn, &x[1], incx); 00141 beta *= rsafmn; 00142 *alpha *= rsafmn; 00143 if (absMACRO(beta) < safmin) { 00144 goto L10; 00145 } 00146 00147 /* New BETA is at most 1, at least SAFMIN */ 00148 00149 i__1 = *n - 1; 00150 xnorm = template_blas_nrm2(&i__1, &x[1], incx); 00151 d__1 = template_lapack_lapy2(alpha, &xnorm); 00152 beta = -template_lapack_d_sign(&d__1, alpha); 00153 *tau = (beta - *alpha) / beta; 00154 i__1 = *n - 1; 00155 d__1 = 1. / (*alpha - beta); 00156 template_blas_scal(&i__1, &d__1, &x[1], incx); 00157 00158 /* If ALPHA is subnormal, it may lose relative accuracy */ 00159 00160 *alpha = beta; 00161 i__1 = knt; 00162 for (j = 1; j <= i__1; ++j) { 00163 *alpha *= safmin; 00164 /* L20: */ 00165 } 00166 } else { 00167 *tau = (beta - *alpha) / beta; 00168 i__1 = *n - 1; 00169 d__1 = 1. / (*alpha - beta); 00170 template_blas_scal(&i__1, &d__1, &x[1], incx); 00171 *alpha = beta; 00172 } 00173 } 00174 00175 return 0; 00176 00177 /* End of DLARFG */ 00178 00179 } /* dlarfg_ */ 00180 00181 #endif