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