ergo
template_lapack_gghrd.h
Go to the documentation of this file.
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