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