00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_LARFG_HEADER 00038 #define TEMPLATE_LAPACK_LARFG_HEADER 00039 00040 #include "template_lapack_lapy2.h" 00041 00042 template<class Treal> 00043 int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, 00044 const integer *incx, Treal *tau) 00045 { 00046 /* -- LAPACK auxiliary routine (version 3.0) -- 00047 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00048 Courant Institute, Argonne National Lab, and Rice University 00049 September 30, 1994 00050 00051 00052 Purpose 00053 ======= 00054 00055 DLARFG generates a real elementary reflector H of order n, such 00056 that 00057 00058 H * ( alpha ) = ( beta ), H' * H = I. 00059 ( x ) ( 0 ) 00060 00061 where alpha and beta are scalars, and x is an (n-1)-element real 00062 vector. H is represented in the form 00063 00064 H = I - tau * ( 1 ) * ( 1 v' ) , 00065 ( v ) 00066 00067 where tau is a real scalar and v is a real (n-1)-element 00068 vector. 00069 00070 If the elements of x are all zero, then tau = 0 and H is taken to be 00071 the unit matrix. 00072 00073 Otherwise 1 <= tau <= 2. 00074 00075 Arguments 00076 ========= 00077 00078 N (input) INTEGER 00079 The order of the elementary reflector. 00080 00081 ALPHA (input/output) DOUBLE PRECISION 00082 On entry, the value alpha. 00083 On exit, it is overwritten with the value beta. 00084 00085 X (input/output) DOUBLE PRECISION array, dimension 00086 (1+(N-2)*abs(INCX)) 00087 On entry, the vector x. 00088 On exit, it is overwritten with the vector v. 00089 00090 INCX (input) INTEGER 00091 The increment between elements of X. INCX > 0. 00092 00093 TAU (output) DOUBLE PRECISION 00094 The value tau. 00095 00096 ===================================================================== 00097 00098 00099 Parameter adjustments */ 00100 /* System generated locals */ 00101 integer i__1; 00102 Treal d__1; 00103 /* Local variables */ 00104 Treal beta; 00105 integer j; 00106 Treal xnorm; 00107 Treal safmin, rsafmn; 00108 integer knt; 00109 00110 --x; 00111 00112 /* Function Body */ 00113 if (*n <= 1) { 00114 *tau = 0.; 00115 return 0; 00116 } 00117 00118 i__1 = *n - 1; 00119 xnorm = template_blas_nrm2(&i__1, &x[1], incx); 00120 00121 if (xnorm == 0.) { 00122 00123 /* H = I */ 00124 00125 *tau = 0.; 00126 } else { 00127 00128 /* general case */ 00129 00130 d__1 = template_lapack_lapy2(alpha, &xnorm); 00131 beta = -template_lapack_d_sign(&d__1, alpha); 00132 safmin = template_lapack_lamch("S", (Treal)0) / template_lapack_lamch("E", (Treal)0); 00133 if (absMACRO(beta) < safmin) { 00134 00135 /* XNORM, BETA may be inaccurate; scale X and recompute them */ 00136 00137 rsafmn = 1. / safmin; 00138 knt = 0; 00139 L10: 00140 ++knt; 00141 i__1 = *n - 1; 00142 template_blas_scal(&i__1, &rsafmn, &x[1], incx); 00143 beta *= rsafmn; 00144 *alpha *= rsafmn; 00145 if (absMACRO(beta) < safmin) { 00146 goto L10; 00147 } 00148 00149 /* New BETA is at most 1, at least SAFMIN */ 00150 00151 i__1 = *n - 1; 00152 xnorm = template_blas_nrm2(&i__1, &x[1], incx); 00153 d__1 = template_lapack_lapy2(alpha, &xnorm); 00154 beta = -template_lapack_d_sign(&d__1, alpha); 00155 *tau = (beta - *alpha) / beta; 00156 i__1 = *n - 1; 00157 d__1 = 1. / (*alpha - beta); 00158 template_blas_scal(&i__1, &d__1, &x[1], incx); 00159 00160 /* If ALPHA is subnormal, it may lose relative accuracy */ 00161 00162 *alpha = beta; 00163 i__1 = knt; 00164 for (j = 1; j <= i__1; ++j) { 00165 *alpha *= safmin; 00166 /* L20: */ 00167 } 00168 } else { 00169 *tau = (beta - *alpha) / beta; 00170 i__1 = *n - 1; 00171 d__1 = 1. / (*alpha - beta); 00172 template_blas_scal(&i__1, &d__1, &x[1], incx); 00173 *alpha = beta; 00174 } 00175 } 00176 00177 return 0; 00178 00179 /* End of DLARFG */ 00180 00181 } /* dlarfg_ */ 00182 00183 #endif