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