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_LARFT_HEADER 00038 #define TEMPLATE_LAPACK_LARFT_HEADER 00039 00040 00041 template<class Treal> 00042 int template_lapack_larft(const char *direct, const char *storev, const integer *n, const integer * 00043 k, Treal *v, const integer *ldv, const Treal *tau, Treal *t, 00044 const integer *ldt) 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 February 29, 1992 00050 00051 00052 Purpose 00053 ======= 00054 00055 DLARFT forms the triangular factor T of a real block reflector H 00056 of order n, which is defined as a product of k elementary reflectors. 00057 00058 If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; 00059 00060 If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. 00061 00062 If STOREV = 'C', the vector which defines the elementary reflector 00063 H(i) is stored in the i-th column of the array V, and 00064 00065 H = I - V * T * V' 00066 00067 If STOREV = 'R', the vector which defines the elementary reflector 00068 H(i) is stored in the i-th row of the array V, and 00069 00070 H = I - V' * T * V 00071 00072 Arguments 00073 ========= 00074 00075 DIRECT (input) CHARACTER*1 00076 Specifies the order in which the elementary reflectors are 00077 multiplied to form the block reflector: 00078 = 'F': H = H(1) H(2) . . . H(k) (Forward) 00079 = 'B': H = H(k) . . . H(2) H(1) (Backward) 00080 00081 STOREV (input) CHARACTER*1 00082 Specifies how the vectors which define the elementary 00083 reflectors are stored (see also Further Details): 00084 = 'C': columnwise 00085 = 'R': rowwise 00086 00087 N (input) INTEGER 00088 The order of the block reflector H. N >= 0. 00089 00090 K (input) INTEGER 00091 The order of the triangular factor T (= the number of 00092 elementary reflectors). K >= 1. 00093 00094 V (input/output) DOUBLE PRECISION array, dimension 00095 (LDV,K) if STOREV = 'C' 00096 (LDV,N) if STOREV = 'R' 00097 The matrix V. See further details. 00098 00099 LDV (input) INTEGER 00100 The leading dimension of the array V. 00101 If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. 00102 00103 TAU (input) DOUBLE PRECISION array, dimension (K) 00104 TAU(i) must contain the scalar factor of the elementary 00105 reflector H(i). 00106 00107 T (output) DOUBLE PRECISION array, dimension (LDT,K) 00108 The k by k triangular factor T of the block reflector. 00109 If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is 00110 lower triangular. The rest of the array is not used. 00111 00112 LDT (input) INTEGER 00113 The leading dimension of the array T. LDT >= K. 00114 00115 Further Details 00116 =============== 00117 00118 The shape of the matrix V and the storage of the vectors which define 00119 the H(i) is best illustrated by the following example with n = 5 and 00120 k = 3. The elements equal to 1 are not stored; the corresponding 00121 array elements are modified but restored on exit. The rest of the 00122 array is not used. 00123 00124 DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': 00125 00126 V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) 00127 ( v1 1 ) ( 1 v2 v2 v2 ) 00128 ( v1 v2 1 ) ( 1 v3 v3 ) 00129 ( v1 v2 v3 ) 00130 ( v1 v2 v3 ) 00131 00132 DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': 00133 00134 V = ( v1 v2 v3 ) V = ( v1 v1 1 ) 00135 ( v1 v2 v3 ) ( v2 v2 v2 1 ) 00136 ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) 00137 ( 1 v3 ) 00138 ( 1 ) 00139 00140 ===================================================================== 00141 00142 00143 Quick return if possible 00144 00145 Parameter adjustments */ 00146 /* Table of constant values */ 00147 integer c__1 = 1; 00148 Treal c_b8 = 0.; 00149 00150 /* System generated locals */ 00151 integer t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3; 00152 Treal d__1; 00153 /* Local variables */ 00154 integer i__, j; 00155 Treal vii; 00156 #define t_ref(a_1,a_2) t[(a_2)*t_dim1 + a_1] 00157 #define v_ref(a_1,a_2) v[(a_2)*v_dim1 + a_1] 00158 00159 00160 v_dim1 = *ldv; 00161 v_offset = 1 + v_dim1 * 1; 00162 v -= v_offset; 00163 --tau; 00164 t_dim1 = *ldt; 00165 t_offset = 1 + t_dim1 * 1; 00166 t -= t_offset; 00167 00168 /* Function Body */ 00169 if (*n == 0) { 00170 return 0; 00171 } 00172 00173 if (template_blas_lsame(direct, "F")) { 00174 i__1 = *k; 00175 for (i__ = 1; i__ <= i__1; ++i__) { 00176 if (tau[i__] == 0.) { 00177 00178 /* H(i) = I */ 00179 00180 i__2 = i__; 00181 for (j = 1; j <= i__2; ++j) { 00182 t_ref(j, i__) = 0.; 00183 /* L10: */ 00184 } 00185 } else { 00186 00187 /* general case */ 00188 00189 vii = v_ref(i__, i__); 00190 v_ref(i__, i__) = 1.; 00191 if (template_blas_lsame(storev, "C")) { 00192 00193 /* T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) */ 00194 00195 i__2 = *n - i__ + 1; 00196 i__3 = i__ - 1; 00197 d__1 = -tau[i__]; 00198 template_blas_gemv("Transpose", &i__2, &i__3, &d__1, &v_ref(i__, 1), 00199 ldv, &v_ref(i__, i__), &c__1, &c_b8, &t_ref(1, 00200 i__), &c__1); 00201 } else { 00202 00203 /* T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' */ 00204 00205 i__2 = i__ - 1; 00206 i__3 = *n - i__ + 1; 00207 d__1 = -tau[i__]; 00208 template_blas_gemv("No transpose", &i__2, &i__3, &d__1, &v_ref(1, i__) 00209 , ldv, &v_ref(i__, i__), ldv, &c_b8, &t_ref(1, 00210 i__), &c__1); 00211 } 00212 v_ref(i__, i__) = vii; 00213 00214 /* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) */ 00215 00216 i__2 = i__ - 1; 00217 template_blas_trmv("Upper", "No transpose", "Non-unit", &i__2, &t[ 00218 t_offset], ldt, &t_ref(1, i__), &c__1); 00219 t_ref(i__, i__) = tau[i__]; 00220 } 00221 /* L20: */ 00222 } 00223 } else { 00224 for (i__ = *k; i__ >= 1; --i__) { 00225 if (tau[i__] == 0.) { 00226 00227 /* H(i) = I */ 00228 00229 i__1 = *k; 00230 for (j = i__; j <= i__1; ++j) { 00231 t_ref(j, i__) = 0.; 00232 /* L30: */ 00233 } 00234 } else { 00235 00236 /* general case */ 00237 00238 if (i__ < *k) { 00239 if (template_blas_lsame(storev, "C")) { 00240 vii = v_ref(*n - *k + i__, i__); 00241 v_ref(*n - *k + i__, i__) = 1.; 00242 00243 /* T(i+1:k,i) := 00244 - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) */ 00245 00246 i__1 = *n - *k + i__; 00247 i__2 = *k - i__; 00248 d__1 = -tau[i__]; 00249 template_blas_gemv("Transpose", &i__1, &i__2, &d__1, &v_ref(1, 00250 i__ + 1), ldv, &v_ref(1, i__), &c__1, &c_b8, & 00251 t_ref(i__ + 1, i__), &c__1); 00252 v_ref(*n - *k + i__, i__) = vii; 00253 } else { 00254 vii = v_ref(i__, *n - *k + i__); 00255 v_ref(i__, *n - *k + i__) = 1.; 00256 00257 /* T(i+1:k,i) := 00258 - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' */ 00259 00260 i__1 = *k - i__; 00261 i__2 = *n - *k + i__; 00262 d__1 = -tau[i__]; 00263 template_blas_gemv("No transpose", &i__1, &i__2, &d__1, &v_ref( 00264 i__ + 1, 1), ldv, &v_ref(i__, 1), ldv, &c_b8, 00265 &t_ref(i__ + 1, i__), &c__1); 00266 v_ref(i__, *n - *k + i__) = vii; 00267 } 00268 00269 /* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) */ 00270 00271 i__1 = *k - i__; 00272 template_blas_trmv("Lower", "No transpose", "Non-unit", &i__1, &t_ref( 00273 i__ + 1, i__ + 1), ldt, &t_ref(i__ + 1, i__), & 00274 c__1); 00275 } 00276 t_ref(i__, i__) = tau[i__]; 00277 } 00278 /* L40: */ 00279 } 00280 } 00281 return 0; 00282 00283 /* End of DLARFT */ 00284 00285 } /* dlarft_ */ 00286 00287 #undef v_ref 00288 #undef t_ref 00289 00290 00291 #endif