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