template_lapack_hgeqz.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_HGEQZ_HEADER
00036 #define TEMPLATE_LAPACK_HGEQZ_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_hgeqz(const char *job, const char *compq, const char *compz, const integer *n, 
00041         const integer *ilo, const integer *ihi, Treal *a, const integer *lda, Treal *
00042         b, const integer *ldb, Treal *alphar, Treal *alphai, Treal *
00043         beta, Treal *q, const integer *ldq, Treal *z__, const integer *ldz, 
00044         Treal *work, const integer *lwork, integer *info)
00045 {
00046 /*  -- LAPACK routine (version 3.0) --   
00047        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00048        Courant Institute, Argonne National Lab, and Rice University   
00049        June 30, 1999   
00050 
00051 
00052     Purpose   
00053     =======   
00054 
00055     DHGEQZ implements a single-/double-shift version of the QZ method for   
00056     finding the generalized eigenvalues   
00057 
00058     w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j)   of the equation   
00059 
00060          det( A - w(i) B ) = 0   
00061 
00062     In addition, the pair A,B may be reduced to generalized Schur form:   
00063     B is upper triangular, and A is block upper triangular, where the   
00064     diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having   
00065     complex generalized eigenvalues (see the description of the argument   
00066     JOB.)   
00067 
00068     If JOB='S', then the pair (A,B) is simultaneously reduced to Schur   
00069     form by applying one orthogonal tranformation (usually called Q) on   
00070     the left and another (usually called Z) on the right.  The 2-by-2   
00071     upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks   
00072     of A will be reduced to positive diagonal matrices.  (I.e.,   
00073     if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and   
00074     B(j+1,j+1) will be positive.)   
00075 
00076     If JOB='E', then at each iteration, the same transformations   
00077     are computed, but they are only applied to those parts of A and B   
00078     which are needed to compute ALPHAR, ALPHAI, and BETAR.   
00079 
00080     If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal   
00081     transformations used to reduce (A,B) are accumulated into the arrays   
00082     Q and Z s.t.:   
00083 
00084          Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*   
00085          Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*   
00086 
00087     Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix   
00088          Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),   
00089          pp. 241--256.   
00090 
00091     Arguments   
00092     =========   
00093 
00094     JOB     (input) CHARACTER*1   
00095             = 'E': compute only ALPHAR, ALPHAI, and BETA.  A and B will   
00096                    not necessarily be put into generalized Schur form.   
00097             = 'S': put A and B into generalized Schur form, as well   
00098                    as computing ALPHAR, ALPHAI, and BETA.   
00099 
00100     COMPQ   (input) CHARACTER*1   
00101             = 'N': do not modify Q.   
00102             = 'V': multiply the array Q on the right by the transpose of   
00103                    the orthogonal tranformation that is applied to the   
00104                    left side of A and B to reduce them to Schur form.   
00105             = 'I': like COMPQ='V', except that Q will be initialized to   
00106                    the identity first.   
00107 
00108     COMPZ   (input) CHARACTER*1   
00109             = 'N': do not modify Z.   
00110             = 'V': multiply the array Z on the right by the orthogonal   
00111                    tranformation that is applied to the right side of   
00112                    A and B to reduce them to Schur form.   
00113             = 'I': like COMPZ='V', except that Z will be initialized to   
00114                    the identity first.   
00115 
00116     N       (input) INTEGER   
00117             The order of the matrices A, B, Q, and Z.  N >= 0.   
00118 
00119     ILO     (input) INTEGER   
00120     IHI     (input) INTEGER   
00121             It is assumed that A is already upper triangular in rows and   
00122             columns 1:ILO-1 and IHI+1:N.   
00123             1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.   
00124 
00125     A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)   
00126             On entry, the N-by-N upper Hessenberg matrix A.  Elements   
00127             below the subdiagonal must be zero.   
00128             If JOB='S', then on exit A and B will have been   
00129                simultaneously reduced to generalized Schur form.   
00130             If JOB='E', then on exit A will have been destroyed.   
00131                The diagonal blocks will be correct, but the off-diagonal   
00132                portion will be meaningless.   
00133 
00134     LDA     (input) INTEGER   
00135             The leading dimension of the array A.  LDA >= max( 1, N ).   
00136 
00137     B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)   
00138             On entry, the N-by-N upper triangular matrix B.  Elements   
00139             below the diagonal must be zero.  2-by-2 blocks in B   
00140             corresponding to 2-by-2 blocks in A will be reduced to   
00141             positive diagonal form.  (I.e., if A(j+1,j) is non-zero,   
00142             then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be   
00143             positive.)   
00144             If JOB='S', then on exit A and B will have been   
00145                simultaneously reduced to Schur form.   
00146             If JOB='E', then on exit B will have been destroyed.   
00147                Elements corresponding to diagonal blocks of A will be   
00148                correct, but the off-diagonal portion will be meaningless.   
00149 
00150     LDB     (input) INTEGER   
00151             The leading dimension of the array B.  LDB >= max( 1, N ).   
00152 
00153     ALPHAR  (output) DOUBLE PRECISION array, dimension (N)   
00154             ALPHAR(1:N) will be set to real parts of the diagonal   
00155             elements of A that would result from reducing A and B to   
00156             Schur form and then further reducing them both to triangular   
00157             form using unitary transformations s.t. the diagonal of B   
00158             was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block   
00159             (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j).   
00160             Note that the (real or complex) values   
00161             (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the   
00162             generalized eigenvalues of the matrix pencil A - wB.   
00163 
00164     ALPHAI  (output) DOUBLE PRECISION array, dimension (N)   
00165             ALPHAI(1:N) will be set to imaginary parts of the diagonal   
00166             elements of A that would result from reducing A and B to   
00167             Schur form and then further reducing them both to triangular   
00168             form using unitary transformations s.t. the diagonal of B   
00169             was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block   
00170             (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0.   
00171             Note that the (real or complex) values   
00172             (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the   
00173             generalized eigenvalues of the matrix pencil A - wB.   
00174 
00175     BETA    (output) DOUBLE PRECISION array, dimension (N)   
00176             BETA(1:N) will be set to the (real) diagonal elements of B   
00177             that would result from reducing A and B to Schur form and   
00178             then further reducing them both to triangular form using   
00179             unitary transformations s.t. the diagonal of B was   
00180             non-negative real.  Thus, if A(j,j) is in a 1-by-1 block   
00181             (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j).   
00182             Note that the (real or complex) values   
00183             (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the   
00184             generalized eigenvalues of the matrix pencil A - wB.   
00185             (Note that BETA(1:N) will always be non-negative, and no   
00186             BETAI is necessary.)   
00187 
00188     Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N)   
00189             If COMPQ='N', then Q will not be referenced.   
00190             If COMPQ='V' or 'I', then the transpose of the orthogonal   
00191                transformations which are applied to A and B on the left   
00192                will be applied to the array Q on the right.   
00193 
00194     LDQ     (input) INTEGER   
00195             The leading dimension of the array Q.  LDQ >= 1.   
00196             If COMPQ='V' or 'I', then LDQ >= N.   
00197 
00198     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)   
00199             If COMPZ='N', then Z will not be referenced.   
00200             If COMPZ='V' or 'I', then the orthogonal transformations   
00201                which are applied to A and B on the right will be applied   
00202                to the array Z on the right.   
00203 
00204     LDZ     (input) INTEGER   
00205             The leading dimension of the array Z.  LDZ >= 1.   
00206             If COMPZ='V' or 'I', then LDZ >= N.   
00207 
00208     WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)   
00209             On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.   
00210 
00211     LWORK   (input) INTEGER   
00212             The dimension of the array WORK.  LWORK >= max(1,N).   
00213 
00214             If LWORK = -1, then a workspace query is assumed; the routine   
00215             only calculates the optimal size of the WORK array, returns   
00216             this value as the first entry of the WORK array, and no error   
00217             message related to LWORK is issued by XERBLA.   
00218 
00219     INFO    (output) INTEGER   
00220             = 0: successful exit   
00221             < 0: if INFO = -i, the i-th argument had an illegal value   
00222             = 1,...,N: the QZ iteration did not converge.  (A,B) is not   
00223                        in Schur form, but ALPHAR(i), ALPHAI(i), and   
00224                        BETA(i), i=INFO+1,...,N should be correct.   
00225             = N+1,...,2*N: the shift calculation failed.  (A,B) is not   
00226                        in Schur form, but ALPHAR(i), ALPHAI(i), and   
00227                        BETA(i), i=INFO-N+1,...,N should be correct.   
00228             > 2*N:     various "impossible" errors.   
00229 
00230     Further Details   
00231     ===============   
00232 
00233     Iteration counters:   
00234 
00235     JITER  -- counts iterations.   
00236     IITER  -- counts iterations run since ILAST was last   
00237               changed.  This is therefore reset only when a 1-by-1 or   
00238               2-by-2 block deflates off the bottom.   
00239 
00240     =====================================================================   
00241 
00242       $                     SAFETY = 1.0E+0 )   
00243 
00244        Decode JOB, COMPQ, COMPZ   
00245 
00246        Parameter adjustments */
00247     /* Table of constant values */
00248      Treal c_b12 = 0.;
00249      Treal c_b13 = 1.;
00250      integer c__1 = 1;
00251      integer c__3 = 3;
00252     
00253     /* System generated locals */
00254     integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 
00255             z_offset, i__1, i__2, i__3, i__4;
00256     Treal d__1, d__2, d__3, d__4;
00257     /* Local variables */
00258      Treal ad11l, ad12l, ad21l, ad22l, ad32l, wabs, atol, btol, 
00259             temp;
00260      Treal temp2, s1inv, c__;
00261      integer j;
00262      Treal s, t, v[3], scale;
00263      integer iiter, ilast, jiter;
00264      Treal anorm, bnorm;
00265      integer maxit;
00266      Treal tempi, tempr, s1, s2, u1, u2;
00267      logical ilazr2;
00268      Treal a11, a12, a21, a22, b11, b22, c12, c21;
00269      integer jc;
00270      Treal an, bn, cl, cq, cr;
00271      integer in;
00272      Treal ascale, bscale, u12, w11;
00273      integer jr;
00274      Treal cz, sl, w12, w21, w22, wi;
00275      Treal sr;
00276      Treal vs, wr;
00277      Treal safmin;
00278      Treal safmax;
00279      Treal eshift;
00280      logical ilschr;
00281      Treal b1a, b2a;
00282      integer icompq, ilastm;
00283      Treal a1i;
00284      integer ischur;
00285      Treal a2i, b1i;
00286      logical ilazro;
00287      integer icompz, ifirst;
00288      Treal b2i;
00289      integer ifrstm;
00290      Treal a1r;
00291      integer istart;
00292      logical ilpivt;
00293      Treal a2r, b1r, b2r;
00294      logical lquery;
00295      Treal wr2, ad11, ad12, ad21, ad22, c11i, c22i;
00296      integer jch;
00297      Treal c11r, c22r, u12l;
00298      logical ilq;
00299      Treal tau, sqi;
00300      logical ilz;
00301      Treal ulp, sqr, szi, szr;
00302 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00303 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00304 #define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]
00305 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
00306 
00307 
00308     a_dim1 = *lda;
00309     a_offset = 1 + a_dim1 * 1;
00310     a -= a_offset;
00311     b_dim1 = *ldb;
00312     b_offset = 1 + b_dim1 * 1;
00313     b -= b_offset;
00314     --alphar;
00315     --alphai;
00316     --beta;
00317     q_dim1 = *ldq;
00318     q_offset = 1 + q_dim1 * 1;
00319     q -= q_offset;
00320     z_dim1 = *ldz;
00321     z_offset = 1 + z_dim1 * 1;
00322     z__ -= z_offset;
00323     --work;
00324 
00325     /* Initialization added by Elias to get rid of compiler warnings. */
00326     ilschr = ilq = ilz = 0;
00327     /* Function Body */
00328     if (template_blas_lsame(job, "E")) {
00329         ilschr = FALSE_;
00330         ischur = 1;
00331     } else if (template_blas_lsame(job, "S")) {
00332         ilschr = TRUE_;
00333         ischur = 2;
00334     } else {
00335         ischur = 0;
00336     }
00337 
00338     if (template_blas_lsame(compq, "N")) {
00339         ilq = FALSE_;
00340         icompq = 1;
00341     } else if (template_blas_lsame(compq, "V")) {
00342         ilq = TRUE_;
00343         icompq = 2;
00344     } else if (template_blas_lsame(compq, "I")) {
00345         ilq = TRUE_;
00346         icompq = 3;
00347     } else {
00348         icompq = 0;
00349     }
00350 
00351     if (template_blas_lsame(compz, "N")) {
00352         ilz = FALSE_;
00353         icompz = 1;
00354     } else if (template_blas_lsame(compz, "V")) {
00355         ilz = TRUE_;
00356         icompz = 2;
00357     } else if (template_blas_lsame(compz, "I")) {
00358         ilz = TRUE_;
00359         icompz = 3;
00360     } else {
00361         icompz = 0;
00362     }
00363 
00364 /*     Check Argument Values */
00365 
00366     *info = 0;
00367     work[1] = (Treal) maxMACRO(1,*n);
00368     lquery = *lwork == -1;
00369     if (ischur == 0) {
00370         *info = -1;
00371     } else if (icompq == 0) {
00372         *info = -2;
00373     } else if (icompz == 0) {
00374         *info = -3;
00375     } else if (*n < 0) {
00376         *info = -4;
00377     } else if (*ilo < 1) {
00378         *info = -5;
00379     } else if (*ihi > *n || *ihi < *ilo - 1) {
00380         *info = -6;
00381     } else if (*lda < *n) {
00382         *info = -8;
00383     } else if (*ldb < *n) {
00384         *info = -10;
00385     } else if (*ldq < 1 || ( ilq && *ldq < *n ) ) {
00386         *info = -15;
00387     } else if (*ldz < 1 || ( ilz && *ldz < *n ) ) {
00388         *info = -17;
00389     } else if (*lwork < maxMACRO(1,*n) && ! lquery) {
00390         *info = -19;
00391     }
00392     if (*info != 0) {
00393         i__1 = -(*info);
00394         template_blas_erbla("HGEQZ ", &i__1);
00395         return 0;
00396     } else if (lquery) {
00397         return 0;
00398     }
00399 
00400 /*     Quick return if possible */
00401 
00402     if (*n <= 0) {
00403         work[1] = 1.;
00404         return 0;
00405     }
00406 
00407 /*     Initialize Q and Z */
00408 
00409     if (icompq == 3) {
00410         template_lapack_laset("Full", n, n, &c_b12, &c_b13, &q[q_offset], ldq);
00411     }
00412     if (icompz == 3) {
00413         template_lapack_laset("Full", n, n, &c_b12, &c_b13, &z__[z_offset], ldz);
00414     }
00415 
00416 /*     Machine Constants */
00417 
00418     in = *ihi + 1 - *ilo;
00419     safmin = template_lapack_lamch("S", (Treal)0);
00420     safmax = 1. / safmin;
00421     ulp = template_lapack_lamch("E", (Treal)0) * template_lapack_lamch("B", (Treal)0);
00422     anorm = dlanhs_("F", &in, &a_ref(*ilo, *ilo), lda, &work[1]);
00423     bnorm = dlanhs_("F", &in, &b_ref(*ilo, *ilo), ldb, &work[1]);
00424 /* Computing MAX */
00425     d__1 = safmin, d__2 = ulp * anorm;
00426     atol = maxMACRO(d__1,d__2);
00427 /* Computing MAX */
00428     d__1 = safmin, d__2 = ulp * bnorm;
00429     btol = maxMACRO(d__1,d__2);
00430     ascale = 1. / maxMACRO(safmin,anorm);
00431     bscale = 1. / maxMACRO(safmin,bnorm);
00432 
00433 /*     Set Eigenvalues IHI+1:N */
00434 
00435     i__1 = *n;
00436     for (j = *ihi + 1; j <= i__1; ++j) {
00437         if (b_ref(j, j) < 0.) {
00438             if (ilschr) {
00439                 i__2 = j;
00440                 for (jr = 1; jr <= i__2; ++jr) {
00441                     a_ref(jr, j) = -a_ref(jr, j);
00442                     b_ref(jr, j) = -b_ref(jr, j);
00443 /* L10: */
00444                 }
00445             } else {
00446                 a_ref(j, j) = -a_ref(j, j);
00447                 b_ref(j, j) = -b_ref(j, j);
00448             }
00449             if (ilz) {
00450                 i__2 = *n;
00451                 for (jr = 1; jr <= i__2; ++jr) {
00452                     z___ref(jr, j) = -z___ref(jr, j);
00453 /* L20: */
00454                 }
00455             }
00456         }
00457         alphar[j] = a_ref(j, j);
00458         alphai[j] = 0.;
00459         beta[j] = b_ref(j, j);
00460 /* L30: */
00461     }
00462 
00463 /*     If IHI < ILO, skip QZ steps */
00464 
00465     if (*ihi < *ilo) {
00466         goto L380;
00467     }
00468 
00469 /*     MAIN QZ ITERATION LOOP   
00470 
00471        Initialize dynamic indices   
00472 
00473        Eigenvalues ILAST+1:N have been found.   
00474           Column operations modify rows IFRSTM:whatever.   
00475           Row operations modify columns whatever:ILASTM.   
00476 
00477        If only eigenvalues are being computed, then   
00478           IFRSTM is the row of the last splitting row above row ILAST;   
00479           this is always at least ILO.   
00480        IITER counts iterations since the last eigenvalue was found,   
00481           to tell when to use an extraordinary shift.   
00482        MAXIT is the maximum number of QZ sweeps allowed. */
00483 
00484     ilast = *ihi;
00485     if (ilschr) {
00486         ifrstm = 1;
00487         ilastm = *n;
00488     } else {
00489         ifrstm = *ilo;
00490         ilastm = *ihi;
00491     }
00492     iiter = 0;
00493     eshift = 0.;
00494     maxit = (*ihi - *ilo + 1) * 30;
00495 
00496     i__1 = maxit;
00497     for (jiter = 1; jiter <= i__1; ++jiter) {
00498 
00499 /*        Split the matrix if possible.   
00500 
00501           Two tests:   
00502              1: A(j,j-1)=0  or  j=ILO   
00503              2: B(j,j)=0 */
00504 
00505         if (ilast == *ilo) {
00506 
00507 /*           Special case: j=ILAST */
00508 
00509             goto L80;
00510         } else {
00511             if ((d__1 = a_ref(ilast, ilast - 1), absMACRO(d__1)) <= atol) {
00512                 a_ref(ilast, ilast - 1) = 0.;
00513                 goto L80;
00514             }
00515         }
00516 
00517         if ((d__1 = b_ref(ilast, ilast), absMACRO(d__1)) <= btol) {
00518             b_ref(ilast, ilast) = 0.;
00519             goto L70;
00520         }
00521 
00522 /*        General case: j<ILAST */
00523 
00524         i__2 = *ilo;
00525         for (j = ilast - 1; j >= i__2; --j) {
00526 
00527 /*           Test 1: for A(j,j-1)=0 or j=ILO */
00528 
00529             if (j == *ilo) {
00530                 ilazro = TRUE_;
00531             } else {
00532                 if ((d__1 = a_ref(j, j - 1), absMACRO(d__1)) <= atol) {
00533                     a_ref(j, j - 1) = 0.;
00534                     ilazro = TRUE_;
00535                 } else {
00536                     ilazro = FALSE_;
00537                 }
00538             }
00539 
00540 /*           Test 2: for B(j,j)=0 */
00541 
00542             if ((d__1 = b_ref(j, j), absMACRO(d__1)) < btol) {
00543                 b_ref(j, j) = 0.;
00544 
00545 /*              Test 1a: Check for 2 consecutive small subdiagonals in A */
00546 
00547                 ilazr2 = FALSE_;
00548                 if (! ilazro) {
00549                     temp = (d__1 = a_ref(j, j - 1), absMACRO(d__1));
00550                     temp2 = (d__1 = a_ref(j, j), absMACRO(d__1));
00551                     tempr = maxMACRO(temp,temp2);
00552                     if (tempr < 1. && tempr != 0.) {
00553                         temp /= tempr;
00554                         temp2 /= tempr;
00555                     }
00556                     if (temp * (ascale * (d__1 = a_ref(j + 1, j), absMACRO(d__1))) 
00557                             <= temp2 * (ascale * atol)) {
00558                         ilazr2 = TRUE_;
00559                     }
00560                 }
00561 
00562 /*              If both tests pass (1 & 2), i.e., the leading diagonal   
00563                 element of B in the block is zero, split a 1x1 block off   
00564                 at the top. (I.e., at the J-th row/column) The leading   
00565                 diagonal element of the remainder can also be zero, so   
00566                 this may have to be done repeatedly. */
00567 
00568                 if (ilazro || ilazr2) {
00569                     i__3 = ilast - 1;
00570                     for (jch = j; jch <= i__3; ++jch) {
00571                         temp = a_ref(jch, jch);
00572                         template_lapack_lartg(&temp, &a_ref(jch + 1, jch), &c__, &s, &a_ref(
00573                                 jch, jch));
00574                         a_ref(jch + 1, jch) = 0.;
00575                         i__4 = ilastm - jch;
00576                         template_blas_rot(&i__4, &a_ref(jch, jch + 1), lda, &a_ref(jch + 
00577                                 1, jch + 1), lda, &c__, &s);
00578                         i__4 = ilastm - jch;
00579                         template_blas_rot(&i__4, &b_ref(jch, jch + 1), ldb, &b_ref(jch + 
00580                                 1, jch + 1), ldb, &c__, &s);
00581                         if (ilq) {
00582                             template_blas_rot(n, &q_ref(1, jch), &c__1, &q_ref(1, jch + 1)
00583                                     , &c__1, &c__, &s);
00584                         }
00585                         if (ilazr2) {
00586                             a_ref(jch, jch - 1) = a_ref(jch, jch - 1) * c__;
00587                         }
00588                         ilazr2 = FALSE_;
00589                         if ((d__1 = b_ref(jch + 1, jch + 1), absMACRO(d__1)) >= 
00590                                 btol) {
00591                             if (jch + 1 >= ilast) {
00592                                 goto L80;
00593                             } else {
00594                                 ifirst = jch + 1;
00595                                 goto L110;
00596                             }
00597                         }
00598                         b_ref(jch + 1, jch + 1) = 0.;
00599 /* L40: */
00600                     }
00601                     goto L70;
00602                 } else {
00603 
00604 /*                 Only test 2 passed -- chase the zero to B(ILAST,ILAST)   
00605                    Then process as in the case B(ILAST,ILAST)=0 */
00606 
00607                     i__3 = ilast - 1;
00608                     for (jch = j; jch <= i__3; ++jch) {
00609                         temp = b_ref(jch, jch + 1);
00610                         template_lapack_lartg(&temp, &b_ref(jch + 1, jch + 1), &c__, &s, &
00611                                 b_ref(jch, jch + 1));
00612                         b_ref(jch + 1, jch + 1) = 0.;
00613                         if (jch < ilastm - 1) {
00614                             i__4 = ilastm - jch - 1;
00615                             template_blas_rot(&i__4, &b_ref(jch, jch + 2), ldb, &b_ref(
00616                                     jch + 1, jch + 2), ldb, &c__, &s);
00617                         }
00618                         i__4 = ilastm - jch + 2;
00619                         template_blas_rot(&i__4, &a_ref(jch, jch - 1), lda, &a_ref(jch + 
00620                                 1, jch - 1), lda, &c__, &s);
00621                         if (ilq) {
00622                             template_blas_rot(n, &q_ref(1, jch), &c__1, &q_ref(1, jch + 1)
00623                                     , &c__1, &c__, &s);
00624                         }
00625                         temp = a_ref(jch + 1, jch);
00626                         template_lapack_lartg(&temp, &a_ref(jch + 1, jch - 1), &c__, &s, &
00627                                 a_ref(jch + 1, jch));
00628                         a_ref(jch + 1, jch - 1) = 0.;
00629                         i__4 = jch + 1 - ifrstm;
00630                         template_blas_rot(&i__4, &a_ref(ifrstm, jch), &c__1, &a_ref(
00631                                 ifrstm, jch - 1), &c__1, &c__, &s);
00632                         i__4 = jch - ifrstm;
00633                         template_blas_rot(&i__4, &b_ref(ifrstm, jch), &c__1, &b_ref(
00634                                 ifrstm, jch - 1), &c__1, &c__, &s);
00635                         if (ilz) {
00636                             template_blas_rot(n, &z___ref(1, jch), &c__1, &z___ref(1, jch 
00637                                     - 1), &c__1, &c__, &s);
00638                         }
00639 /* L50: */
00640                     }
00641                     goto L70;
00642                 }
00643             } else if (ilazro) {
00644 
00645 /*              Only test 1 passed -- work on J:ILAST */
00646 
00647                 ifirst = j;
00648                 goto L110;
00649             }
00650 
00651 /*           Neither test passed -- try next J   
00652 
00653    L60: */
00654         }
00655 
00656 /*        (Drop-through is "impossible") */
00657 
00658         *info = *n + 1;
00659         goto L420;
00660 
00661 /*        B(ILAST,ILAST)=0 -- clear A(ILAST,ILAST-1) to split off a   
00662           1x1 block. */
00663 
00664 L70:
00665         temp = a_ref(ilast, ilast);
00666         template_lapack_lartg(&temp, &a_ref(ilast, ilast - 1), &c__, &s, &a_ref(ilast, 
00667                 ilast));
00668         a_ref(ilast, ilast - 1) = 0.;
00669         i__2 = ilast - ifrstm;
00670         template_blas_rot(&i__2, &a_ref(ifrstm, ilast), &c__1, &a_ref(ifrstm, ilast - 1), 
00671                 &c__1, &c__, &s);
00672         i__2 = ilast - ifrstm;
00673         template_blas_rot(&i__2, &b_ref(ifrstm, ilast), &c__1, &b_ref(ifrstm, ilast - 1), 
00674                 &c__1, &c__, &s);
00675         if (ilz) {
00676             template_blas_rot(n, &z___ref(1, ilast), &c__1, &z___ref(1, ilast - 1), &c__1,
00677                      &c__, &s);
00678         }
00679 
00680 /*        A(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,   
00681                                 and BETA */
00682 
00683 L80:
00684         if (b_ref(ilast, ilast) < 0.) {
00685             if (ilschr) {
00686                 i__2 = ilast;
00687                 for (j = ifrstm; j <= i__2; ++j) {
00688                     a_ref(j, ilast) = -a_ref(j, ilast);
00689                     b_ref(j, ilast) = -b_ref(j, ilast);
00690 /* L90: */
00691                 }
00692             } else {
00693                 a_ref(ilast, ilast) = -a_ref(ilast, ilast);
00694                 b_ref(ilast, ilast) = -b_ref(ilast, ilast);
00695             }
00696             if (ilz) {
00697                 i__2 = *n;
00698                 for (j = 1; j <= i__2; ++j) {
00699                     z___ref(j, ilast) = -z___ref(j, ilast);
00700 /* L100: */
00701                 }
00702             }
00703         }
00704         alphar[ilast] = a_ref(ilast, ilast);
00705         alphai[ilast] = 0.;
00706         beta[ilast] = b_ref(ilast, ilast);
00707 
00708 /*        Go to next block -- exit if finished. */
00709 
00710         --ilast;
00711         if (ilast < *ilo) {
00712             goto L380;
00713         }
00714 
00715 /*        Reset counters */
00716 
00717         iiter = 0;
00718         eshift = 0.;
00719         if (! ilschr) {
00720             ilastm = ilast;
00721             if (ifrstm > ilast) {
00722                 ifrstm = *ilo;
00723             }
00724         }
00725         goto L350;
00726 
00727 /*        QZ step   
00728 
00729           This iteration only involves rows/columns IFIRST:ILAST. We   
00730           assume IFIRST < ILAST, and that the diagonal of B is non-zero. */
00731 
00732 L110:
00733         ++iiter;
00734         if (! ilschr) {
00735             ifrstm = ifirst;
00736         }
00737 
00738 /*        Compute single shifts.   
00739 
00740           At this point, IFIRST < ILAST, and the diagonal elements of   
00741           B(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in   
00742           magnitude) */
00743 
00744         if (iiter / 10 * 10 == iiter) {
00745 
00746 /*           Exceptional shift.  Chosen for no particularly good reason.   
00747              (Single shift only.) */
00748 
00749             if ((Treal) maxit * safmin * (d__1 = a_ref(ilast - 1, ilast),
00750                      absMACRO(d__1)) < (d__2 = b_ref(ilast - 1, ilast - 1), absMACRO(
00751                     d__2))) {
00752                 eshift += a_ref(ilast - 1, ilast) / b_ref(ilast - 1, ilast - 
00753                         1);
00754             } else {
00755                 eshift += 1. / (safmin * (Treal) maxit);
00756             }
00757             s1 = 1.;
00758             wr = eshift;
00759 
00760         } else {
00761 
00762 /*           Shifts based on the generalized eigenvalues of the   
00763              bottom-right 2x2 block of A and B. The first eigenvalue   
00764              returned by DLAG2 is the Wilkinson shift (AEP p.512), */
00765 
00766             d__1 = safmin * 100.;
00767             template_lapack_lag2(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast 
00768                     - 1), ldb, &d__1, &s1, &s2, &wr, &wr2, &wi);
00769 
00770 /* Computing MAX   
00771    Computing MAX */
00772             d__3 = 1., d__4 = absMACRO(wr), d__3 = maxMACRO(d__3,d__4), d__4 = absMACRO(wi);
00773             d__1 = s1, d__2 = safmin * maxMACRO(d__3,d__4);
00774             temp = maxMACRO(d__1,d__2);
00775             if (wi != 0.) {
00776                 goto L200;
00777             }
00778         }
00779 
00780 /*        Fiddle with shift to avoid overflow */
00781 
00782         temp = minMACRO(ascale,1.) * (safmax * .5);
00783         if (s1 > temp) {
00784             scale = temp / s1;
00785         } else {
00786             scale = 1.;
00787         }
00788 
00789         temp = minMACRO(bscale,1.) * (safmax * .5);
00790         if (absMACRO(wr) > temp) {
00791 /* Computing MIN */
00792             d__1 = scale, d__2 = temp / absMACRO(wr);
00793             scale = minMACRO(d__1,d__2);
00794         }
00795         s1 = scale * s1;
00796         wr = scale * wr;
00797 
00798 /*        Now check for two consecutive small subdiagonals. */
00799 
00800         i__2 = ifirst + 1;
00801         for (j = ilast - 1; j >= i__2; --j) {
00802             istart = j;
00803             temp = (d__1 = s1 * a_ref(j, j - 1), absMACRO(d__1));
00804             temp2 = (d__1 = s1 * a_ref(j, j) - wr * b_ref(j, j), absMACRO(d__1));
00805             tempr = maxMACRO(temp,temp2);
00806             if (tempr < 1. && tempr != 0.) {
00807                 temp /= tempr;
00808                 temp2 /= tempr;
00809             }
00810             if ((d__1 = ascale * a_ref(j + 1, j) * temp, absMACRO(d__1)) <= ascale 
00811                     * atol * temp2) {
00812                 goto L130;
00813             }
00814 /* L120: */
00815         }
00816 
00817         istart = ifirst;
00818 L130:
00819 
00820 /*        Do an implicit single-shift QZ sweep.   
00821 
00822           Initial Q */
00823 
00824         temp = s1 * a_ref(istart, istart) - wr * b_ref(istart, istart);
00825         temp2 = s1 * a_ref(istart + 1, istart);
00826         template_lapack_lartg(&temp, &temp2, &c__, &s, &tempr);
00827 
00828 /*        Sweep */
00829 
00830         i__2 = ilast - 1;
00831         for (j = istart; j <= i__2; ++j) {
00832             if (j > istart) {
00833                 temp = a_ref(j, j - 1);
00834                 template_lapack_lartg(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j - 
00835                         1));
00836                 a_ref(j + 1, j - 1) = 0.;
00837             }
00838 
00839             i__3 = ilastm;
00840             for (jc = j; jc <= i__3; ++jc) {
00841                 temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc);
00842                 a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc);
00843                 a_ref(j, jc) = temp;
00844                 temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc);
00845                 b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc);
00846                 b_ref(j, jc) = temp2;
00847 /* L140: */
00848             }
00849             if (ilq) {
00850                 i__3 = *n;
00851                 for (jr = 1; jr <= i__3; ++jr) {
00852                     temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1);
00853                     q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j 
00854                             + 1);
00855                     q_ref(jr, j) = temp;
00856 /* L150: */
00857                 }
00858             }
00859 
00860             temp = b_ref(j + 1, j + 1);
00861             template_lapack_lartg(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1));
00862             b_ref(j + 1, j) = 0.;
00863 
00864 /* Computing MIN */
00865             i__4 = j + 2;
00866             i__3 = minMACRO(i__4,ilast);
00867             for (jr = ifrstm; jr <= i__3; ++jr) {
00868                 temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j);
00869                 a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j);
00870                 a_ref(jr, j + 1) = temp;
00871 /* L160: */
00872             }
00873             i__3 = j;
00874             for (jr = ifrstm; jr <= i__3; ++jr) {
00875                 temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j);
00876                 b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j);
00877                 b_ref(jr, j + 1) = temp;
00878 /* L170: */
00879             }
00880             if (ilz) {
00881                 i__3 = *n;
00882                 for (jr = 1; jr <= i__3; ++jr) {
00883                     temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j);
00884                     z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref(
00885                             jr, j);
00886                     z___ref(jr, j + 1) = temp;
00887 /* L180: */
00888                 }
00889             }
00890 /* L190: */
00891         }
00892 
00893         goto L350;
00894 
00895 /*        Use Francis double-shift   
00896 
00897           Note: the Francis double-shift should work with real shifts,   
00898                 but only if the block is at least 3x3.   
00899                 This code may break if this point is reached with   
00900                 a 2x2 block with real eigenvalues. */
00901 
00902 L200:
00903         if (ifirst + 1 == ilast) {
00904 
00905 /*           Special case -- 2x2 block with complex eigenvectors   
00906 
00907              Step 1: Standardize, that is, rotate so that   
00908 
00909                          ( B11  0  )   
00910                      B = (         )  with B11 non-negative.   
00911                          (  0  B22 ) */
00912 
00913             template_lapack_lasv2(&b_ref(ilast - 1, ilast - 1), &b_ref(ilast - 1, ilast), &
00914                     b_ref(ilast, ilast), &b22, &b11, &sr, &cr, &sl, &cl);
00915 
00916             if (b11 < 0.) {
00917                 cr = -cr;
00918                 sr = -sr;
00919                 b11 = -b11;
00920                 b22 = -b22;
00921             }
00922 
00923             i__2 = ilastm + 1 - ifirst;
00924             template_blas_rot(&i__2, &a_ref(ilast - 1, ilast - 1), lda, &a_ref(ilast, 
00925                     ilast - 1), lda, &cl, &sl);
00926             i__2 = ilast + 1 - ifrstm;
00927             template_blas_rot(&i__2, &a_ref(ifrstm, ilast - 1), &c__1, &a_ref(ifrstm, 
00928                     ilast), &c__1, &cr, &sr);
00929 
00930             if (ilast < ilastm) {
00931                 i__2 = ilastm - ilast;
00932                 template_blas_rot(&i__2, &b_ref(ilast - 1, ilast + 1), ldb, &b_ref(ilast, 
00933                         ilast + 1), lda, &cl, &sl);
00934             }
00935             if (ifrstm < ilast - 1) {
00936                 i__2 = ifirst - ifrstm;
00937                 template_blas_rot(&i__2, &b_ref(ifrstm, ilast - 1), &c__1, &b_ref(ifrstm, 
00938                         ilast), &c__1, &cr, &sr);
00939             }
00940 
00941             if (ilq) {
00942                 template_blas_rot(n, &q_ref(1, ilast - 1), &c__1, &q_ref(1, ilast), &c__1,
00943                          &cl, &sl);
00944             }
00945             if (ilz) {
00946                 template_blas_rot(n, &z___ref(1, ilast - 1), &c__1, &z___ref(1, ilast), &
00947                         c__1, &cr, &sr);
00948             }
00949 
00950             b_ref(ilast - 1, ilast - 1) = b11;
00951             b_ref(ilast - 1, ilast) = 0.;
00952             b_ref(ilast, ilast - 1) = 0.;
00953             b_ref(ilast, ilast) = b22;
00954 
00955 /*           If B22 is negative, negate column ILAST */
00956 
00957             if (b22 < 0.) {
00958                 i__2 = ilast;
00959                 for (j = ifrstm; j <= i__2; ++j) {
00960                     a_ref(j, ilast) = -a_ref(j, ilast);
00961                     b_ref(j, ilast) = -b_ref(j, ilast);
00962 /* L210: */
00963                 }
00964 
00965                 if (ilz) {
00966                     i__2 = *n;
00967                     for (j = 1; j <= i__2; ++j) {
00968                         z___ref(j, ilast) = -z___ref(j, ilast);
00969 /* L220: */
00970                     }
00971                 }
00972             }
00973 
00974 /*           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)   
00975 
00976              Recompute shift */
00977 
00978             d__1 = safmin * 100.;
00979             template_lapack_lag2(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast 
00980                     - 1), ldb, &d__1, &s1, &temp, &wr, &temp2, &wi);
00981 
00982 /*           If standardization has perturbed the shift onto real line,   
00983              do another (real single-shift) QR step. */
00984 
00985             if (wi == 0.) {
00986                 goto L350;
00987             }
00988             s1inv = 1. / s1;
00989 
00990 /*           Do EISPACK (QZVAL) computation of alpha and beta */
00991 
00992             a11 = a_ref(ilast - 1, ilast - 1);
00993             a21 = a_ref(ilast, ilast - 1);
00994             a12 = a_ref(ilast - 1, ilast);
00995             a22 = a_ref(ilast, ilast);
00996 
00997 /*           Compute complex Givens rotation on right   
00998              (Assume some element of C = (sA - wB) > unfl )   
00999                               __   
01000              (sA - wB) ( CZ   -SZ )   
01001                        ( SZ    CZ ) */
01002 
01003             c11r = s1 * a11 - wr * b11;
01004             c11i = -wi * b11;
01005             c12 = s1 * a12;
01006             c21 = s1 * a21;
01007             c22r = s1 * a22 - wr * b22;
01008             c22i = -wi * b22;
01009 
01010             if (absMACRO(c11r) + absMACRO(c11i) + absMACRO(c12) > absMACRO(c21) + absMACRO(c22r) + absMACRO(
01011                     c22i)) {
01012                 t = template_lapack_lapy3(&c12, &c11r, &c11i);
01013                 cz = c12 / t;
01014                 szr = -c11r / t;
01015                 szi = -c11i / t;
01016             } else {
01017                 cz = template_lapack_lapy2(&c22r, &c22i);
01018                 if (cz <= safmin) {
01019                     cz = 0.;
01020                     szr = 1.;
01021                     szi = 0.;
01022                 } else {
01023                     tempr = c22r / cz;
01024                     tempi = c22i / cz;
01025                     t = template_lapack_lapy2(&cz, &c21);
01026                     cz /= t;
01027                     szr = -c21 * tempr / t;
01028                     szi = c21 * tempi / t;
01029                 }
01030             }
01031 
01032 /*           Compute Givens rotation on left   
01033 
01034              (  CQ   SQ )   
01035              (  __      )  A or B   
01036              ( -SQ   CQ ) */
01037 
01038             an = absMACRO(a11) + absMACRO(a12) + absMACRO(a21) + absMACRO(a22);
01039             bn = absMACRO(b11) + absMACRO(b22);
01040             wabs = absMACRO(wr) + absMACRO(wi);
01041             if (s1 * an > wabs * bn) {
01042                 cq = cz * b11;
01043                 sqr = szr * b22;
01044                 sqi = -szi * b22;
01045             } else {
01046                 a1r = cz * a11 + szr * a12;
01047                 a1i = szi * a12;
01048                 a2r = cz * a21 + szr * a22;
01049                 a2i = szi * a22;
01050                 cq = template_lapack_lapy2(&a1r, &a1i);
01051                 if (cq <= safmin) {
01052                     cq = 0.;
01053                     sqr = 1.;
01054                     sqi = 0.;
01055                 } else {
01056                     tempr = a1r / cq;
01057                     tempi = a1i / cq;
01058                     sqr = tempr * a2r + tempi * a2i;
01059                     sqi = tempi * a2r - tempr * a2i;
01060                 }
01061             }
01062             t = template_lapack_lapy3(&cq, &sqr, &sqi);
01063             cq /= t;
01064             sqr /= t;
01065             sqi /= t;
01066 
01067 /*           Compute diagonal elements of QBZ */
01068 
01069             tempr = sqr * szr - sqi * szi;
01070             tempi = sqr * szi + sqi * szr;
01071             b1r = cq * cz * b11 + tempr * b22;
01072             b1i = tempi * b22;
01073             b1a = template_lapack_lapy2(&b1r, &b1i);
01074             b2r = cq * cz * b22 + tempr * b11;
01075             b2i = -tempi * b11;
01076             b2a = template_lapack_lapy2(&b2r, &b2i);
01077 
01078 /*           Normalize so beta > 0, and Im( alpha1 ) > 0 */
01079 
01080             beta[ilast - 1] = b1a;
01081             beta[ilast] = b2a;
01082             alphar[ilast - 1] = wr * b1a * s1inv;
01083             alphai[ilast - 1] = wi * b1a * s1inv;
01084             alphar[ilast] = wr * b2a * s1inv;
01085             alphai[ilast] = -(wi * b2a) * s1inv;
01086 
01087 /*           Step 3: Go to next block -- exit if finished. */
01088 
01089             ilast = ifirst - 1;
01090             if (ilast < *ilo) {
01091                 goto L380;
01092             }
01093 
01094 /*           Reset counters */
01095 
01096             iiter = 0;
01097             eshift = 0.;
01098             if (! ilschr) {
01099                 ilastm = ilast;
01100                 if (ifrstm > ilast) {
01101                     ifrstm = *ilo;
01102                 }
01103             }
01104             goto L350;
01105         } else {
01106 
01107 /*           Usual case: 3x3 or larger block, using Francis implicit   
01108                          double-shift   
01109 
01110                                       2   
01111              Eigenvalue equation is  w  - c w + d = 0,   
01112 
01113                                            -1 2        -1   
01114              so compute 1st column of  (A B  )  - c A B   + d   
01115              using the formula in QZIT (from EISPACK)   
01116 
01117              We assume that the block is at least 3x3 */
01118 
01119             ad11 = ascale * a_ref(ilast - 1, ilast - 1) / (bscale * b_ref(
01120                     ilast - 1, ilast - 1));
01121             ad21 = ascale * a_ref(ilast, ilast - 1) / (bscale * b_ref(ilast - 
01122                     1, ilast - 1));
01123             ad12 = ascale * a_ref(ilast - 1, ilast) / (bscale * b_ref(ilast, 
01124                     ilast));
01125             ad22 = ascale * a_ref(ilast, ilast) / (bscale * b_ref(ilast, 
01126                     ilast));
01127             u12 = b_ref(ilast - 1, ilast) / b_ref(ilast, ilast);
01128             ad11l = ascale * a_ref(ifirst, ifirst) / (bscale * b_ref(ifirst, 
01129                     ifirst));
01130             ad21l = ascale * a_ref(ifirst + 1, ifirst) / (bscale * b_ref(
01131                     ifirst, ifirst));
01132             ad12l = ascale * a_ref(ifirst, ifirst + 1) / (bscale * b_ref(
01133                     ifirst + 1, ifirst + 1));
01134             ad22l = ascale * a_ref(ifirst + 1, ifirst + 1) / (bscale * b_ref(
01135                     ifirst + 1, ifirst + 1));
01136             ad32l = ascale * a_ref(ifirst + 2, ifirst + 1) / (bscale * b_ref(
01137                     ifirst + 1, ifirst + 1));
01138             u12l = b_ref(ifirst, ifirst + 1) / b_ref(ifirst + 1, ifirst + 1);
01139 
01140             v[0] = (ad11 - ad11l) * (ad22 - ad11l) - ad12 * ad21 + ad21 * u12 
01141                     * ad11l + (ad12l - ad11l * u12l) * ad21l;
01142             v[1] = (ad22l - ad11l - ad21l * u12l - (ad11 - ad11l) - (ad22 - 
01143                     ad11l) + ad21 * u12) * ad21l;
01144             v[2] = ad32l * ad21l;
01145 
01146             istart = ifirst;
01147 
01148             template_lapack_larfg(&c__3, v, &v[1], &c__1, &tau);
01149             v[0] = 1.;
01150 
01151 /*           Sweep */
01152 
01153             i__2 = ilast - 2;
01154             for (j = istart; j <= i__2; ++j) {
01155 
01156 /*              All but last elements: use 3x3 Householder transforms.   
01157 
01158                 Zero (j-1)st column of A */
01159 
01160                 if (j > istart) {
01161                     v[0] = a_ref(j, j - 1);
01162                     v[1] = a_ref(j + 1, j - 1);
01163                     v[2] = a_ref(j + 2, j - 1);
01164 
01165                     template_lapack_larfg(&c__3, &a_ref(j, j - 1), &v[1], &c__1, &tau);
01166                     v[0] = 1.;
01167                     a_ref(j + 1, j - 1) = 0.;
01168                     a_ref(j + 2, j - 1) = 0.;
01169                 }
01170 
01171                 i__3 = ilastm;
01172                 for (jc = j; jc <= i__3; ++jc) {
01173                     temp = tau * (a_ref(j, jc) + v[1] * a_ref(j + 1, jc) + v[
01174                             2] * a_ref(j + 2, jc));
01175                     a_ref(j, jc) = a_ref(j, jc) - temp;
01176                     a_ref(j + 1, jc) = a_ref(j + 1, jc) - temp * v[1];
01177                     a_ref(j + 2, jc) = a_ref(j + 2, jc) - temp * v[2];
01178                     temp2 = tau * (b_ref(j, jc) + v[1] * b_ref(j + 1, jc) + v[
01179                             2] * b_ref(j + 2, jc));
01180                     b_ref(j, jc) = b_ref(j, jc) - temp2;
01181                     b_ref(j + 1, jc) = b_ref(j + 1, jc) - temp2 * v[1];
01182                     b_ref(j + 2, jc) = b_ref(j + 2, jc) - temp2 * v[2];
01183 /* L230: */
01184                 }
01185                 if (ilq) {
01186                     i__3 = *n;
01187                     for (jr = 1; jr <= i__3; ++jr) {
01188                         temp = tau * (q_ref(jr, j) + v[1] * q_ref(jr, j + 1) 
01189                                 + v[2] * q_ref(jr, j + 2));
01190                         q_ref(jr, j) = q_ref(jr, j) - temp;
01191                         q_ref(jr, j + 1) = q_ref(jr, j + 1) - temp * v[1];
01192                         q_ref(jr, j + 2) = q_ref(jr, j + 2) - temp * v[2];
01193 /* L240: */
01194                     }
01195                 }
01196 
01197 /*              Zero j-th column of B (see DLAGBC for details)   
01198 
01199                 Swap rows to pivot */
01200 
01201                 ilpivt = FALSE_;
01202 /* Computing MAX */
01203                 d__3 = (d__1 = b_ref(j + 1, j + 1), absMACRO(d__1)), d__4 = (d__2 =
01204                          b_ref(j + 1, j + 2), absMACRO(d__2));
01205                 temp = maxMACRO(d__3,d__4);
01206 /* Computing MAX */
01207                 d__3 = (d__1 = b_ref(j + 2, j + 1), absMACRO(d__1)), d__4 = (d__2 =
01208                          b_ref(j + 2, j + 2), absMACRO(d__2));
01209                 temp2 = maxMACRO(d__3,d__4);
01210                 if (maxMACRO(temp,temp2) < safmin) {
01211                     scale = 0.;
01212                     u1 = 1.;
01213                     u2 = 0.;
01214                     goto L250;
01215                 } else if (temp >= temp2) {
01216                     w11 = b_ref(j + 1, j + 1);
01217                     w21 = b_ref(j + 2, j + 1);
01218                     w12 = b_ref(j + 1, j + 2);
01219                     w22 = b_ref(j + 2, j + 2);
01220                     u1 = b_ref(j + 1, j);
01221                     u2 = b_ref(j + 2, j);
01222                 } else {
01223                     w21 = b_ref(j + 1, j + 1);
01224                     w11 = b_ref(j + 2, j + 1);
01225                     w22 = b_ref(j + 1, j + 2);
01226                     w12 = b_ref(j + 2, j + 2);
01227                     u2 = b_ref(j + 1, j);
01228                     u1 = b_ref(j + 2, j);
01229                 }
01230 
01231 /*              Swap columns if nec. */
01232 
01233                 if (absMACRO(w12) > absMACRO(w11)) {
01234                     ilpivt = TRUE_;
01235                     temp = w12;
01236                     temp2 = w22;
01237                     w12 = w11;
01238                     w22 = w21;
01239                     w11 = temp;
01240                     w21 = temp2;
01241                 }
01242 
01243 /*              LU-factor */
01244 
01245                 temp = w21 / w11;
01246                 u2 -= temp * u1;
01247                 w22 -= temp * w12;
01248                 w21 = 0.;
01249 
01250 /*              Compute SCALE */
01251 
01252                 scale = 1.;
01253                 if (absMACRO(w22) < safmin) {
01254                     scale = 0.;
01255                     u2 = 1.;
01256                     u1 = -w12 / w11;
01257                     goto L250;
01258                 }
01259                 if (absMACRO(w22) < absMACRO(u2)) {
01260                     scale = (d__1 = w22 / u2, absMACRO(d__1));
01261                 }
01262                 if (absMACRO(w11) < absMACRO(u1)) {
01263 /* Computing MIN */
01264                     d__2 = scale, d__3 = (d__1 = w11 / u1, absMACRO(d__1));
01265                     scale = minMACRO(d__2,d__3);
01266                 }
01267 
01268 /*              Solve */
01269 
01270                 u2 = scale * u2 / w22;
01271                 u1 = (scale * u1 - w12 * u2) / w11;
01272 
01273 L250:
01274                 if (ilpivt) {
01275                     temp = u2;
01276                     u2 = u1;
01277                     u1 = temp;
01278                 }
01279 
01280 /*              Compute Householder Vector   
01281 
01282    Computing 2nd power */
01283                 d__1 = scale;
01284 /* Computing 2nd power */
01285                 d__2 = u1;
01286 /* Computing 2nd power */
01287                 d__3 = u2;
01288                 t = template_blas_sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
01289                 tau = scale / t + 1.;
01290                 vs = -1. / (scale + t);
01291                 v[0] = 1.;
01292                 v[1] = vs * u1;
01293                 v[2] = vs * u2;
01294 
01295 /*              Apply transformations from the right.   
01296 
01297    Computing MIN */
01298                 i__4 = j + 3;
01299                 i__3 = minMACRO(i__4,ilast);
01300                 for (jr = ifrstm; jr <= i__3; ++jr) {
01301                     temp = tau * (a_ref(jr, j) + v[1] * a_ref(jr, j + 1) + v[
01302                             2] * a_ref(jr, j + 2));
01303                     a_ref(jr, j) = a_ref(jr, j) - temp;
01304                     a_ref(jr, j + 1) = a_ref(jr, j + 1) - temp * v[1];
01305                     a_ref(jr, j + 2) = a_ref(jr, j + 2) - temp * v[2];
01306 /* L260: */
01307                 }
01308                 i__3 = j + 2;
01309                 for (jr = ifrstm; jr <= i__3; ++jr) {
01310                     temp = tau * (b_ref(jr, j) + v[1] * b_ref(jr, j + 1) + v[
01311                             2] * b_ref(jr, j + 2));
01312                     b_ref(jr, j) = b_ref(jr, j) - temp;
01313                     b_ref(jr, j + 1) = b_ref(jr, j + 1) - temp * v[1];
01314                     b_ref(jr, j + 2) = b_ref(jr, j + 2) - temp * v[2];
01315 /* L270: */
01316                 }
01317                 if (ilz) {
01318                     i__3 = *n;
01319                     for (jr = 1; jr <= i__3; ++jr) {
01320                         temp = tau * (z___ref(jr, j) + v[1] * z___ref(jr, j + 
01321                                 1) + v[2] * z___ref(jr, j + 2));
01322                         z___ref(jr, j) = z___ref(jr, j) - temp;
01323                         z___ref(jr, j + 1) = z___ref(jr, j + 1) - temp * v[1];
01324                         z___ref(jr, j + 2) = z___ref(jr, j + 2) - temp * v[2];
01325 /* L280: */
01326                     }
01327                 }
01328                 b_ref(j + 1, j) = 0.;
01329                 b_ref(j + 2, j) = 0.;
01330 /* L290: */
01331             }
01332 
01333 /*           Last elements: Use Givens rotations   
01334 
01335              Rotations from the left */
01336 
01337             j = ilast - 1;
01338             temp = a_ref(j, j - 1);
01339             template_lapack_lartg(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j - 1));
01340             a_ref(j + 1, j - 1) = 0.;
01341 
01342             i__2 = ilastm;
01343             for (jc = j; jc <= i__2; ++jc) {
01344                 temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc);
01345                 a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc);
01346                 a_ref(j, jc) = temp;
01347                 temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc);
01348                 b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc);
01349                 b_ref(j, jc) = temp2;
01350 /* L300: */
01351             }
01352             if (ilq) {
01353                 i__2 = *n;
01354                 for (jr = 1; jr <= i__2; ++jr) {
01355                     temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1);
01356                     q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j 
01357                             + 1);
01358                     q_ref(jr, j) = temp;
01359 /* L310: */
01360                 }
01361             }
01362 
01363 /*           Rotations from the right. */
01364 
01365             temp = b_ref(j + 1, j + 1);
01366             template_lapack_lartg(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1));
01367             b_ref(j + 1, j) = 0.;
01368 
01369             i__2 = ilast;
01370             for (jr = ifrstm; jr <= i__2; ++jr) {
01371                 temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j);
01372                 a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j);
01373                 a_ref(jr, j + 1) = temp;
01374 /* L320: */
01375             }
01376             i__2 = ilast - 1;
01377             for (jr = ifrstm; jr <= i__2; ++jr) {
01378                 temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j);
01379                 b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j);
01380                 b_ref(jr, j + 1) = temp;
01381 /* L330: */
01382             }
01383             if (ilz) {
01384                 i__2 = *n;
01385                 for (jr = 1; jr <= i__2; ++jr) {
01386                     temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j);
01387                     z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref(
01388                             jr, j);
01389                     z___ref(jr, j + 1) = temp;
01390 /* L340: */
01391                 }
01392             }
01393 
01394 /*           End of Double-Shift code */
01395 
01396         }
01397 
01398         goto L350;
01399 
01400 /*        End of iteration loop */
01401 
01402 L350:
01403 /* L360: */
01404         ;
01405     }
01406 
01407 /*     Drop-through = non-convergence   
01408 
01409    L370: */
01410     *info = ilast;
01411     goto L420;
01412 
01413 /*     Successful completion of all QZ steps */
01414 
01415 L380:
01416 
01417 /*     Set Eigenvalues 1:ILO-1 */
01418 
01419     i__1 = *ilo - 1;
01420     for (j = 1; j <= i__1; ++j) {
01421         if (b_ref(j, j) < 0.) {
01422             if (ilschr) {
01423                 i__2 = j;
01424                 for (jr = 1; jr <= i__2; ++jr) {
01425                     a_ref(jr, j) = -a_ref(jr, j);
01426                     b_ref(jr, j) = -b_ref(jr, j);
01427 /* L390: */
01428                 }
01429             } else {
01430                 a_ref(j, j) = -a_ref(j, j);
01431                 b_ref(j, j) = -b_ref(j, j);
01432             }
01433             if (ilz) {
01434                 i__2 = *n;
01435                 for (jr = 1; jr <= i__2; ++jr) {
01436                     z___ref(jr, j) = -z___ref(jr, j);
01437 /* L400: */
01438                 }
01439             }
01440         }
01441         alphar[j] = a_ref(j, j);
01442         alphai[j] = 0.;
01443         beta[j] = b_ref(j, j);
01444 /* L410: */
01445     }
01446 
01447 /*     Normal Termination */
01448 
01449     *info = 0;
01450 
01451 /*     Exit (other than argument error) -- return optimal workspace size */
01452 
01453 L420:
01454     work[1] = (Treal) (*n);
01455     return 0;
01456 
01457 /*     End of DHGEQZ */
01458 
01459 } /* dhgeqz_ */
01460 
01461 #undef z___ref
01462 #undef q_ref
01463 #undef b_ref
01464 #undef a_ref
01465 
01466 
01467 #endif

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1