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_ORGTR_HEADER 00038 #define TEMPLATE_LAPACK_ORGTR_HEADER 00039 00040 00041 template<class Treal> 00042 int template_lapack_orgtr(const char *uplo, const integer *n, Treal *a, const integer * 00043 lda, const 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 DORGTR generates a real orthogonal matrix Q which is defined as the 00055 product of n-1 elementary reflectors of order N, as returned by 00056 DSYTRD: 00057 00058 if UPLO = 'U', Q = H(n-1) . . . H(2) H(1), 00059 00060 if UPLO = 'L', Q = H(1) H(2) . . . H(n-1). 00061 00062 Arguments 00063 ========= 00064 00065 UPLO (input) CHARACTER*1 00066 = 'U': Upper triangle of A contains elementary reflectors 00067 from DSYTRD; 00068 = 'L': Lower triangle of A contains elementary reflectors 00069 from DSYTRD. 00070 00071 N (input) INTEGER 00072 The order of the matrix Q. N >= 0. 00073 00074 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00075 On entry, the vectors which define the elementary reflectors, 00076 as returned by DSYTRD. 00077 On exit, the N-by-N orthogonal matrix Q. 00078 00079 LDA (input) INTEGER 00080 The leading dimension of the array A. LDA >= max(1,N). 00081 00082 TAU (input) DOUBLE PRECISION array, dimension (N-1) 00083 TAU(i) must contain the scalar factor of the elementary 00084 reflector H(i), as returned by DSYTRD. 00085 00086 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 00087 On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00088 00089 LWORK (input) INTEGER 00090 The dimension of the array WORK. LWORK >= max(1,N-1). 00091 For optimum performance LWORK >= (N-1)*NB, where NB is 00092 the optimal blocksize. 00093 00094 If LWORK = -1, then a workspace query is assumed; the routine 00095 only calculates the optimal size of the WORK array, returns 00096 this value as the first entry of the WORK array, and no error 00097 message related to LWORK is issued by XERBLA. 00098 00099 INFO (output) INTEGER 00100 = 0: successful exit 00101 < 0: if INFO = -i, the i-th argument had an illegal value 00102 00103 ===================================================================== 00104 00105 00106 Test the input arguments 00107 00108 Parameter adjustments */ 00109 /* Table of constant values */ 00110 integer c__1 = 1; 00111 integer c_n1 = -1; 00112 00113 /* System generated locals */ 00114 integer a_dim1, a_offset, i__1, i__2, i__3; 00115 /* Local variables */ 00116 integer i__, j; 00117 integer iinfo; 00118 logical upper; 00119 integer nb; 00120 integer lwkopt; 00121 logical lquery; 00122 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00123 00124 00125 a_dim1 = *lda; 00126 a_offset = 1 + a_dim1 * 1; 00127 a -= a_offset; 00128 --tau; 00129 --work; 00130 00131 /* Initialization added by Elias to get rid of compiler warnings. */ 00132 lwkopt = 0; 00133 /* Function Body */ 00134 *info = 0; 00135 lquery = *lwork == -1; 00136 upper = template_blas_lsame(uplo, "U"); 00137 if (! upper && ! template_blas_lsame(uplo, "L")) { 00138 *info = -1; 00139 } else if (*n < 0) { 00140 *info = -2; 00141 } else if (*lda < maxMACRO(1,*n)) { 00142 *info = -4; 00143 } else /* if(complicated condition) */ { 00144 /* Computing MAX */ 00145 i__1 = 1, i__2 = *n - 1; 00146 if (*lwork < maxMACRO(i__1,i__2) && ! lquery) { 00147 *info = -7; 00148 } 00149 } 00150 00151 if (*info == 0) { 00152 if (upper) { 00153 i__1 = *n - 1; 00154 i__2 = *n - 1; 00155 i__3 = *n - 1; 00156 nb = template_lapack_ilaenv(&c__1, "DORGQL", " ", &i__1, &i__2, &i__3, &c_n1, ( 00157 ftnlen)6, (ftnlen)1); 00158 } else { 00159 i__1 = *n - 1; 00160 i__2 = *n - 1; 00161 i__3 = *n - 1; 00162 nb = template_lapack_ilaenv(&c__1, "DORGQR", " ", &i__1, &i__2, &i__3, &c_n1, ( 00163 ftnlen)6, (ftnlen)1); 00164 } 00165 /* Computing MAX */ 00166 i__1 = 1, i__2 = *n - 1; 00167 lwkopt = maxMACRO(i__1,i__2) * nb; 00168 work[1] = (Treal) lwkopt; 00169 } 00170 00171 if (*info != 0) { 00172 i__1 = -(*info); 00173 template_blas_erbla("ORGTR ", &i__1); 00174 return 0; 00175 } else if (lquery) { 00176 return 0; 00177 } 00178 00179 /* Quick return if possible */ 00180 00181 if (*n == 0) { 00182 work[1] = 1.; 00183 return 0; 00184 } 00185 00186 if (upper) { 00187 00188 /* Q was determined by a call to DSYTRD with UPLO = 'U' 00189 00190 Shift the vectors which define the elementary reflectors one 00191 column to the left, and set the last row and column of Q to 00192 those of the unit matrix */ 00193 00194 i__1 = *n - 1; 00195 for (j = 1; j <= i__1; ++j) { 00196 i__2 = j - 1; 00197 for (i__ = 1; i__ <= i__2; ++i__) { 00198 a_ref(i__, j) = a_ref(i__, j + 1); 00199 /* L10: */ 00200 } 00201 a_ref(*n, j) = 0.; 00202 /* L20: */ 00203 } 00204 i__1 = *n - 1; 00205 for (i__ = 1; i__ <= i__1; ++i__) { 00206 a_ref(i__, *n) = 0.; 00207 /* L30: */ 00208 } 00209 a_ref(*n, *n) = 1.; 00210 00211 /* Generate Q(1:n-1,1:n-1) */ 00212 00213 i__1 = *n - 1; 00214 i__2 = *n - 1; 00215 i__3 = *n - 1; 00216 template_lapack_orgql(&i__1, &i__2, &i__3, &a[a_offset], lda, &tau[1], &work[1], 00217 lwork, &iinfo); 00218 00219 } else { 00220 00221 /* Q was determined by a call to DSYTRD with UPLO = 'L'. 00222 00223 Shift the vectors which define the elementary reflectors one 00224 column to the right, and set the first row and column of Q to 00225 those of the unit matrix */ 00226 00227 for (j = *n; j >= 2; --j) { 00228 a_ref(1, j) = 0.; 00229 i__1 = *n; 00230 for (i__ = j + 1; i__ <= i__1; ++i__) { 00231 a_ref(i__, j) = a_ref(i__, j - 1); 00232 /* L40: */ 00233 } 00234 /* L50: */ 00235 } 00236 a_ref(1, 1) = 1.; 00237 i__1 = *n; 00238 for (i__ = 2; i__ <= i__1; ++i__) { 00239 a_ref(i__, 1) = 0.; 00240 /* L60: */ 00241 } 00242 if (*n > 1) { 00243 00244 /* Generate Q(2:n,2:n) */ 00245 00246 i__1 = *n - 1; 00247 i__2 = *n - 1; 00248 i__3 = *n - 1; 00249 template_lapack_orgqr( 00250 &i__1, 00251 &i__2, 00252 &i__3, 00253 &a_ref(2, 2), 00254 lda, 00255 &tau[1], 00256 &work[1], 00257 lwork, 00258 &iinfo 00259 ); 00260 } 00261 } 00262 work[1] = (Treal) lwkopt; 00263 return 0; 00264 00265 /* End of DORGTR */ 00266 00267 } /* dorgtr_ */ 00268 00269 #undef a_ref 00270 00271 00272 #endif