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_GEQRF_HEADER 00038 #define TEMPLATE_LAPACK_GEQRF_HEADER 00039 00040 00041 template<class Treal> 00042 int template_lapack_geqrf(const integer *m, const integer *n, Treal *a, const integer * 00043 lda, Treal *tau, Treal *work, const integer *lwork, integer *info) 00044 { 00045 /* -- LAPACK routine (version 3.0) -- 00046 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00047 Courant Institute, Argonne National Lab, and Rice University 00048 June 30, 1999 00049 00050 00051 Purpose 00052 ======= 00053 00054 DGEQRF computes a QR factorization of a real M-by-N matrix A: 00055 A = Q * R. 00056 00057 Arguments 00058 ========= 00059 00060 M (input) INTEGER 00061 The number of rows of the matrix A. M >= 0. 00062 00063 N (input) INTEGER 00064 The number of columns of the matrix A. N >= 0. 00065 00066 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00067 On entry, the M-by-N matrix A. 00068 On exit, the elements on and above the diagonal of the array 00069 contain the min(M,N)-by-N upper trapezoidal matrix R (R is 00070 upper triangular if m >= n); the elements below the diagonal, 00071 with the array TAU, represent the orthogonal matrix Q as a 00072 product of min(m,n) elementary reflectors (see Further 00073 Details). 00074 00075 LDA (input) INTEGER 00076 The leading dimension of the array A. LDA >= max(1,M). 00077 00078 TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) 00079 The scalar factors of the elementary reflectors (see Further 00080 Details). 00081 00082 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 00083 On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00084 00085 LWORK (input) INTEGER 00086 The dimension of the array WORK. LWORK >= max(1,N). 00087 For optimum performance LWORK >= N*NB, where NB is 00088 the optimal blocksize. 00089 00090 If LWORK = -1, then a workspace query is assumed; the routine 00091 only calculates the optimal size of the WORK array, returns 00092 this value as the first entry of the WORK array, and no error 00093 message related to LWORK is issued by XERBLA. 00094 00095 INFO (output) INTEGER 00096 = 0: successful exit 00097 < 0: if INFO = -i, the i-th argument had an illegal value 00098 00099 Further Details 00100 =============== 00101 00102 The matrix Q is represented as a product of elementary reflectors 00103 00104 Q = H(1) H(2) . . . H(k), where k = min(m,n). 00105 00106 Each H(i) has the form 00107 00108 H(i) = I - tau * v * v' 00109 00110 where tau is a real scalar, and v is a real vector with 00111 v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), 00112 and tau in TAU(i). 00113 00114 ===================================================================== 00115 00116 00117 Test the input arguments 00118 00119 Parameter adjustments */ 00120 /* Table of constant values */ 00121 integer c__1 = 1; 00122 integer c_n1 = -1; 00123 integer c__3 = 3; 00124 integer c__2 = 2; 00125 00126 /* System generated locals */ 00127 integer a_dim1, a_offset, i__1, i__2, i__3, i__4; 00128 /* Local variables */ 00129 integer i__, k, nbmin, iinfo; 00130 integer ib, nb; 00131 integer nx; 00132 integer ldwork, lwkopt; 00133 logical lquery; 00134 integer iws; 00135 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00136 00137 00138 a_dim1 = *lda; 00139 a_offset = 1 + a_dim1 * 1; 00140 a -= a_offset; 00141 --tau; 00142 --work; 00143 00144 /* Function Body */ 00145 *info = 0; 00146 nb = template_lapack_ilaenv(&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen) 00147 1); 00148 lwkopt = *n * nb; 00149 work[1] = (Treal) lwkopt; 00150 lquery = *lwork == -1; 00151 if (*m < 0) { 00152 *info = -1; 00153 } else if (*n < 0) { 00154 *info = -2; 00155 } else if (*lda < maxMACRO(1,*m)) { 00156 *info = -4; 00157 } else if (*lwork < maxMACRO(1,*n) && ! lquery) { 00158 *info = -7; 00159 } 00160 if (*info != 0) { 00161 i__1 = -(*info); 00162 template_blas_erbla("GEQRF ", &i__1); 00163 return 0; 00164 } else if (lquery) { 00165 return 0; 00166 } 00167 00168 /* Quick return if possible */ 00169 00170 k = minMACRO(*m,*n); 00171 if (k == 0) { 00172 work[1] = 1.; 00173 return 0; 00174 } 00175 00176 nbmin = 2; 00177 nx = 0; 00178 iws = *n; 00179 if (nb > 1 && nb < k) { 00180 00181 /* Determine when to cross over from blocked to unblocked code. 00182 00183 Computing MAX */ 00184 i__1 = 0, i__2 = template_lapack_ilaenv(&c__3, "DGEQRF", " ", m, n, &c_n1, &c_n1, ( 00185 ftnlen)6, (ftnlen)1); 00186 nx = maxMACRO(i__1,i__2); 00187 if (nx < k) { 00188 00189 /* Determine if workspace is large enough for blocked code. */ 00190 00191 ldwork = *n; 00192 iws = ldwork * nb; 00193 if (*lwork < iws) { 00194 00195 /* Not enough workspace to use optimal NB: reduce NB and 00196 determine the minimum value of NB. */ 00197 00198 nb = *lwork / ldwork; 00199 /* Computing MAX */ 00200 i__1 = 2, i__2 = template_lapack_ilaenv(&c__2, "DGEQRF", " ", m, n, &c_n1, & 00201 c_n1, (ftnlen)6, (ftnlen)1); 00202 nbmin = maxMACRO(i__1,i__2); 00203 } 00204 } 00205 } 00206 00207 if (nb >= nbmin && nb < k && nx < k) { 00208 00209 /* Use blocked code initially */ 00210 00211 i__1 = k - nx; 00212 i__2 = nb; 00213 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) { 00214 /* Computing MIN */ 00215 i__3 = k - i__ + 1; 00216 ib = minMACRO(i__3,nb); 00217 00218 /* Compute the QR factorization of the current block 00219 A(i:m,i:i+ib-1) */ 00220 00221 i__3 = *m - i__ + 1; 00222 template_lapack_geqr2(&i__3, &ib, &a_ref(i__, i__), lda, &tau[i__], &work[1], & 00223 iinfo); 00224 if (i__ + ib <= *n) { 00225 00226 /* Form the triangular factor of the block reflector 00227 H = H(i) H(i+1) . . . H(i+ib-1) */ 00228 00229 i__3 = *m - i__ + 1; 00230 template_lapack_larft("Forward", "Columnwise", &i__3, &ib, &a_ref(i__, i__), 00231 lda, &tau[i__], &work[1], &ldwork); 00232 00233 /* Apply H' to A(i:m,i+ib:n) from the left */ 00234 00235 i__3 = *m - i__ + 1; 00236 i__4 = *n - i__ - ib + 1; 00237 template_lapack_larfb("Left", "Transpose", "Forward", "Columnwise", &i__3, & 00238 i__4, &ib, &a_ref(i__, i__), lda, &work[1], &ldwork, & 00239 a_ref(i__, i__ + ib), lda, &work[ib + 1], &ldwork); 00240 } 00241 /* L10: */ 00242 } 00243 } else { 00244 i__ = 1; 00245 } 00246 00247 /* Use unblocked code to factor the last or only block. */ 00248 00249 if (i__ <= k) { 00250 i__2 = *m - i__ + 1; 00251 i__1 = *n - i__ + 1; 00252 template_lapack_geqr2(&i__2, &i__1, &a_ref(i__, i__), lda, &tau[i__], &work[1], & 00253 iinfo); 00254 } 00255 00256 work[1] = (Treal) iws; 00257 return 0; 00258 00259 /* End of DGEQRF */ 00260 00261 } /* dgeqrf_ */ 00262 00263 #undef a_ref 00264 00265 00266 #endif