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_ORM2R_HEADER 00036 #define TEMPLATE_LAPACK_ORM2R_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_orm2r(const char *side, const char *trans, const integer *m, const integer *n, 00041 const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal * 00042 c__, const integer *ldc, Treal *work, integer *info) 00043 { 00044 /* -- LAPACK 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 DORM2R overwrites the general real m by n matrix C with 00054 00055 Q * C if SIDE = 'L' and TRANS = 'N', or 00056 00057 Q'* C if SIDE = 'L' and TRANS = 'T', or 00058 00059 C * Q if SIDE = 'R' and TRANS = 'N', or 00060 00061 C * Q' if SIDE = 'R' and TRANS = 'T', 00062 00063 where Q is a real orthogonal matrix defined as the product of k 00064 elementary reflectors 00065 00066 Q = H(1) H(2) . . . H(k) 00067 00068 as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n 00069 if SIDE = 'R'. 00070 00071 Arguments 00072 ========= 00073 00074 SIDE (input) CHARACTER*1 00075 = 'L': apply Q or Q' from the Left 00076 = 'R': apply Q or Q' from the Right 00077 00078 TRANS (input) CHARACTER*1 00079 = 'N': apply Q (No transpose) 00080 = 'T': apply Q' (Transpose) 00081 00082 M (input) INTEGER 00083 The number of rows of the matrix C. M >= 0. 00084 00085 N (input) INTEGER 00086 The number of columns of the matrix C. N >= 0. 00087 00088 K (input) INTEGER 00089 The number of elementary reflectors whose product defines 00090 the matrix Q. 00091 If SIDE = 'L', M >= K >= 0; 00092 if SIDE = 'R', N >= K >= 0. 00093 00094 A (input) DOUBLE PRECISION array, dimension (LDA,K) 00095 The i-th column must contain the vector which defines the 00096 elementary reflector H(i), for i = 1,2,...,k, as returned by 00097 DGEQRF in the first k columns of its array argument A. 00098 A is modified by the routine but restored on exit. 00099 00100 LDA (input) INTEGER 00101 The leading dimension of the array A. 00102 If SIDE = 'L', LDA >= max(1,M); 00103 if SIDE = 'R', LDA >= max(1,N). 00104 00105 TAU (input) DOUBLE PRECISION array, dimension (K) 00106 TAU(i) must contain the scalar factor of the elementary 00107 reflector H(i), as returned by DGEQRF. 00108 00109 C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 00110 On entry, the m by n matrix C. 00111 On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. 00112 00113 LDC (input) INTEGER 00114 The leading dimension of the array C. LDC >= max(1,M). 00115 00116 WORK (workspace) DOUBLE PRECISION array, dimension 00117 (N) if SIDE = 'L', 00118 (M) if SIDE = 'R' 00119 00120 INFO (output) INTEGER 00121 = 0: successful exit 00122 < 0: if INFO = -i, the i-th argument had an illegal value 00123 00124 ===================================================================== 00125 00126 00127 Test the input arguments 00128 00129 Parameter adjustments */ 00130 /* Table of constant values */ 00131 integer c__1 = 1; 00132 00133 /* System generated locals */ 00134 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2; 00135 /* Local variables */ 00136 logical left; 00137 integer i__; 00138 integer i1, i2, i3, ic, jc, mi, ni, nq; 00139 logical notran; 00140 Treal aii; 00141 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00142 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1] 00143 00144 00145 a_dim1 = *lda; 00146 a_offset = 1 + a_dim1 * 1; 00147 a -= a_offset; 00148 --tau; 00149 c_dim1 = *ldc; 00150 c_offset = 1 + c_dim1 * 1; 00151 c__ -= c_offset; 00152 --work; 00153 00154 /* Function Body */ 00155 *info = 0; 00156 left = template_blas_lsame(side, "L"); 00157 notran = template_blas_lsame(trans, "N"); 00158 00159 /* NQ is the order of Q */ 00160 00161 if (left) { 00162 nq = *m; 00163 } else { 00164 nq = *n; 00165 } 00166 if (! left && ! template_blas_lsame(side, "R")) { 00167 *info = -1; 00168 } else if (! notran && ! template_blas_lsame(trans, "T")) { 00169 *info = -2; 00170 } else if (*m < 0) { 00171 *info = -3; 00172 } else if (*n < 0) { 00173 *info = -4; 00174 } else if (*k < 0 || *k > nq) { 00175 *info = -5; 00176 } else if (*lda < maxMACRO(1,nq)) { 00177 *info = -7; 00178 } else if (*ldc < maxMACRO(1,*m)) { 00179 *info = -10; 00180 } 00181 if (*info != 0) { 00182 i__1 = -(*info); 00183 template_blas_erbla("ORM2R ", &i__1); 00184 return 0; 00185 } 00186 00187 /* Quick return if possible */ 00188 00189 if (*m == 0 || *n == 0 || *k == 0) { 00190 return 0; 00191 } 00192 00193 if ( ( left && ! notran ) || ( ! left && notran ) ) { 00194 i1 = 1; 00195 i2 = *k; 00196 i3 = 1; 00197 } else { 00198 i1 = *k; 00199 i2 = 1; 00200 i3 = -1; 00201 } 00202 00203 if (left) { 00204 ni = *n; 00205 jc = 1; 00206 } else { 00207 mi = *m; 00208 ic = 1; 00209 } 00210 00211 i__1 = i2; 00212 i__2 = i3; 00213 for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) { 00214 if (left) { 00215 00216 /* H(i) is applied to C(i:m,1:n) */ 00217 00218 mi = *m - i__ + 1; 00219 ic = i__; 00220 } else { 00221 00222 /* H(i) is applied to C(1:m,i:n) */ 00223 00224 ni = *n - i__ + 1; 00225 jc = i__; 00226 } 00227 00228 /* Apply H(i) */ 00229 00230 aii = a_ref(i__, i__); 00231 a_ref(i__, i__) = 1.; 00232 template_lapack_larf(side, &mi, &ni, &a_ref(i__, i__), &c__1, &tau[i__], &c___ref( 00233 ic, jc), ldc, &work[1]); 00234 a_ref(i__, i__) = aii; 00235 /* L10: */ 00236 } 00237 return 0; 00238 00239 /* End of DORM2R */ 00240 00241 } /* dorm2r_ */ 00242 00243 #undef c___ref 00244 #undef a_ref 00245 00246 00247 #endif