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_GGHRD_HEADER 00036 #define TEMPLATE_LAPACK_GGHRD_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_gghrd(const char *compq, const char *compz, const integer *n, const integer * 00041 ilo, const integer *ihi, Treal *a, const integer *lda, Treal *b, 00042 const integer *ldb, Treal *q, const integer *ldq, Treal *z__, const integer * 00043 ldz, 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 September 30, 1994 00049 00050 00051 Purpose 00052 ======= 00053 00054 DGGHRD reduces a pair of real matrices (A,B) to generalized upper 00055 Hessenberg form using orthogonal transformations, where A is a 00056 general matrix and B is upper triangular: Q' * A * Z = H and 00057 Q' * B * Z = T, where H is upper Hessenberg, T is upper triangular, 00058 and Q and Z are orthogonal, and ' means transpose. 00059 00060 The orthogonal matrices Q and Z are determined as products of Givens 00061 rotations. They may either be formed explicitly, or they may be 00062 postmultiplied into input matrices Q1 and Z1, so that 00063 00064 Q1 * A * Z1' = (Q1*Q) * H * (Z1*Z)' 00065 Q1 * B * Z1' = (Q1*Q) * T * (Z1*Z)' 00066 00067 Arguments 00068 ========= 00069 00070 COMPQ (input) CHARACTER*1 00071 = 'N': do not compute Q; 00072 = 'I': Q is initialized to the unit matrix, and the 00073 orthogonal matrix Q is returned; 00074 = 'V': Q must contain an orthogonal matrix Q1 on entry, 00075 and the product Q1*Q is returned. 00076 00077 COMPZ (input) CHARACTER*1 00078 = 'N': do not compute Z; 00079 = 'I': Z is initialized to the unit matrix, and the 00080 orthogonal matrix Z is returned; 00081 = 'V': Z must contain an orthogonal matrix Z1 on entry, 00082 and the product Z1*Z is returned. 00083 00084 N (input) INTEGER 00085 The order of the matrices A and B. N >= 0. 00086 00087 ILO (input) INTEGER 00088 IHI (input) INTEGER 00089 It is assumed that A is already upper triangular in rows and 00090 columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally set 00091 by a previous call to DGGBAL; otherwise they should be set 00092 to 1 and N respectively. 00093 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 00094 00095 A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00096 On entry, the N-by-N general matrix to be reduced. 00097 On exit, the upper triangle and the first subdiagonal of A 00098 are overwritten with the upper Hessenberg matrix H, and the 00099 rest is set to zero. 00100 00101 LDA (input) INTEGER 00102 The leading dimension of the array A. LDA >= max(1,N). 00103 00104 B (input/output) DOUBLE PRECISION array, dimension (LDB, N) 00105 On entry, the N-by-N upper triangular matrix B. 00106 On exit, the upper triangular matrix T = Q' B Z. The 00107 elements below the diagonal are set to zero. 00108 00109 LDB (input) INTEGER 00110 The leading dimension of the array B. LDB >= max(1,N). 00111 00112 Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) 00113 If COMPQ='N': Q is not referenced. 00114 If COMPQ='I': on entry, Q need not be set, and on exit it 00115 contains the orthogonal matrix Q, where Q' 00116 is the product of the Givens transformations 00117 which are applied to A and B on the left. 00118 If COMPQ='V': on entry, Q must contain an orthogonal matrix 00119 Q1, and on exit this is overwritten by Q1*Q. 00120 00121 LDQ (input) INTEGER 00122 The leading dimension of the array Q. 00123 LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. 00124 00125 Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) 00126 If COMPZ='N': Z is not referenced. 00127 If COMPZ='I': on entry, Z need not be set, and on exit it 00128 contains the orthogonal matrix Z, which is 00129 the product of the Givens transformations 00130 which are applied to A and B on the right. 00131 If COMPZ='V': on entry, Z must contain an orthogonal matrix 00132 Z1, and on exit this is overwritten by Z1*Z. 00133 00134 LDZ (input) INTEGER 00135 The leading dimension of the array Z. 00136 LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. 00137 00138 INFO (output) INTEGER 00139 = 0: successful exit. 00140 < 0: if INFO = -i, the i-th argument had an illegal value. 00141 00142 Further Details 00143 =============== 00144 00145 This routine reduces A to Hessenberg and B to triangular form by 00146 an unblocked reduction, as described in _Matrix_Computations_, 00147 by Golub and Van Loan (Johns Hopkins Press.) 00148 00149 ===================================================================== 00150 00151 00152 Decode COMPQ 00153 00154 Parameter adjustments */ 00155 /* Table of constant values */ 00156 Treal c_b10 = 0.; 00157 Treal c_b11 = 1.; 00158 integer c__1 = 1; 00159 00160 /* System generated locals */ 00161 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 00162 z_offset, i__1, i__2, i__3; 00163 /* Local variables */ 00164 integer jcol; 00165 Treal temp; 00166 integer jrow; 00167 Treal c__, s; 00168 integer icompq, icompz; 00169 logical ilq, ilz; 00170 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00171 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] 00172 #define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1] 00173 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1] 00174 00175 00176 a_dim1 = *lda; 00177 a_offset = 1 + a_dim1 * 1; 00178 a -= a_offset; 00179 b_dim1 = *ldb; 00180 b_offset = 1 + b_dim1 * 1; 00181 b -= b_offset; 00182 q_dim1 = *ldq; 00183 q_offset = 1 + q_dim1 * 1; 00184 q -= q_offset; 00185 z_dim1 = *ldz; 00186 z_offset = 1 + z_dim1 * 1; 00187 z__ -= z_offset; 00188 00189 /* Initialization added by Elias to get rid of compiler warnings. */ 00190 ilq = ilz = 0; 00191 /* Function Body */ 00192 if (template_blas_lsame(compq, "N")) { 00193 ilq = FALSE_; 00194 icompq = 1; 00195 } else if (template_blas_lsame(compq, "V")) { 00196 ilq = TRUE_; 00197 icompq = 2; 00198 } else if (template_blas_lsame(compq, "I")) { 00199 ilq = TRUE_; 00200 icompq = 3; 00201 } else { 00202 icompq = 0; 00203 } 00204 00205 /* Decode COMPZ */ 00206 00207 if (template_blas_lsame(compz, "N")) { 00208 ilz = FALSE_; 00209 icompz = 1; 00210 } else if (template_blas_lsame(compz, "V")) { 00211 ilz = TRUE_; 00212 icompz = 2; 00213 } else if (template_blas_lsame(compz, "I")) { 00214 ilz = TRUE_; 00215 icompz = 3; 00216 } else { 00217 icompz = 0; 00218 } 00219 00220 /* Test the input parameters. */ 00221 00222 *info = 0; 00223 if (icompq <= 0) { 00224 *info = -1; 00225 } else if (icompz <= 0) { 00226 *info = -2; 00227 } else if (*n < 0) { 00228 *info = -3; 00229 } else if (*ilo < 1) { 00230 *info = -4; 00231 } else if (*ihi > *n || *ihi < *ilo - 1) { 00232 *info = -5; 00233 } else if (*lda < maxMACRO(1,*n)) { 00234 *info = -7; 00235 } else if (*ldb < maxMACRO(1,*n)) { 00236 *info = -9; 00237 } else if ( ( ilq && *ldq < *n ) || *ldq < 1) { 00238 *info = -11; 00239 } else if ( ( ilz && *ldz < *n ) || *ldz < 1) { 00240 *info = -13; 00241 } 00242 if (*info != 0) { 00243 i__1 = -(*info); 00244 template_blas_erbla("GGHRD ", &i__1); 00245 return 0; 00246 } 00247 00248 /* Initialize Q and Z if desired. */ 00249 00250 if (icompq == 3) { 00251 template_lapack_laset("Full", n, n, &c_b10, &c_b11, &q[q_offset], ldq); 00252 } 00253 if (icompz == 3) { 00254 template_lapack_laset("Full", n, n, &c_b10, &c_b11, &z__[z_offset], ldz); 00255 } 00256 00257 /* Quick return if possible */ 00258 00259 if (*n <= 1) { 00260 return 0; 00261 } 00262 00263 /* Zero out lower triangle of B */ 00264 00265 i__1 = *n - 1; 00266 for (jcol = 1; jcol <= i__1; ++jcol) { 00267 i__2 = *n; 00268 for (jrow = jcol + 1; jrow <= i__2; ++jrow) { 00269 b_ref(jrow, jcol) = 0.; 00270 /* L10: */ 00271 } 00272 /* L20: */ 00273 } 00274 00275 /* Reduce A and B */ 00276 00277 i__1 = *ihi - 2; 00278 for (jcol = *ilo; jcol <= i__1; ++jcol) { 00279 00280 i__2 = jcol + 2; 00281 for (jrow = *ihi; jrow >= i__2; --jrow) { 00282 00283 /* Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) */ 00284 00285 temp = a_ref(jrow - 1, jcol); 00286 template_lapack_lartg(&temp, &a_ref(jrow, jcol), &c__, &s, &a_ref(jrow - 1, 00287 jcol)); 00288 a_ref(jrow, jcol) = 0.; 00289 i__3 = *n - jcol; 00290 template_blas_rot(&i__3, &a_ref(jrow - 1, jcol + 1), lda, &a_ref(jrow, jcol + 00291 1), lda, &c__, &s); 00292 i__3 = *n + 2 - jrow; 00293 template_blas_rot(&i__3, &b_ref(jrow - 1, jrow - 1), ldb, &b_ref(jrow, jrow - 00294 1), ldb, &c__, &s); 00295 if (ilq) { 00296 template_blas_rot(n, &q_ref(1, jrow - 1), &c__1, &q_ref(1, jrow), &c__1, & 00297 c__, &s); 00298 } 00299 00300 /* Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) */ 00301 00302 temp = b_ref(jrow, jrow); 00303 template_lapack_lartg(&temp, &b_ref(jrow, jrow - 1), &c__, &s, &b_ref(jrow, 00304 jrow)); 00305 b_ref(jrow, jrow - 1) = 0.; 00306 template_blas_rot(ihi, &a_ref(1, jrow), &c__1, &a_ref(1, jrow - 1), &c__1, & 00307 c__, &s); 00308 i__3 = jrow - 1; 00309 template_blas_rot(&i__3, &b_ref(1, jrow), &c__1, &b_ref(1, jrow - 1), &c__1, & 00310 c__, &s); 00311 if (ilz) { 00312 template_blas_rot(n, &z___ref(1, jrow), &c__1, &z___ref(1, jrow - 1), & 00313 c__1, &c__, &s); 00314 } 00315 /* L30: */ 00316 } 00317 /* L40: */ 00318 } 00319 00320 return 0; 00321 00322 /* End of DGGHRD */ 00323 00324 } /* dgghrd_ */ 00325 00326 #undef z___ref 00327 #undef q_ref 00328 #undef b_ref 00329 #undef a_ref 00330 00331 00332 #endif