template_lapack_tgevc.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_TGEVC_HEADER
00036 #define TEMPLATE_LAPACK_TGEVC_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_tgevc(const char *side, const char *howmny, const logical *select, 
00041         const integer *n, const Treal *a, const integer *lda, const Treal *b, const integer *ldb, 
00042         Treal *vl, const integer *ldvl, Treal *vr, const integer *ldvr, const integer 
00043         *mm, integer *m, Treal *work, 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        June 30, 1999   
00049 
00050 
00051 
00052     Purpose   
00053     =======   
00054 
00055     DTGEVC computes some or all of the right and/or left generalized   
00056     eigenvectors of a pair of real upper triangular matrices (A,B).   
00057 
00058     The right generalized eigenvector x and the left generalized   
00059     eigenvector y of (A,B) corresponding to a generalized eigenvalue   
00060     w are defined by:   
00061 
00062             (A - wB) * x = 0  and  y**H * (A - wB) = 0   
00063 
00064     where y**H denotes the conjugate tranpose of y.   
00065 
00066     If an eigenvalue w is determined by zero diagonal elements of both A   
00067     and B, a unit vector is returned as the corresponding eigenvector.   
00068 
00069     If all eigenvectors are requested, the routine may either return   
00070     the matrices X and/or Y of right or left eigenvectors of (A,B), or   
00071     the products Z*X and/or Q*Y, where Z and Q are input orthogonal   
00072     matrices.  If (A,B) was obtained from the generalized real-Schur   
00073     factorization of an original pair of matrices   
00074        (A0,B0) = (Q*A*Z**H,Q*B*Z**H),   
00075     then Z*X and Q*Y are the matrices of right or left eigenvectors of   
00076     A.   
00077 
00078     A must be block upper triangular, with 1-by-1 and 2-by-2 diagonal   
00079     blocks.  Corresponding to each 2-by-2 diagonal block is a complex   
00080     conjugate pair of eigenvalues and eigenvectors; only one   
00081     eigenvector of the pair is computed, namely the one corresponding   
00082     to the eigenvalue with positive imaginary part.   
00083 
00084     Arguments   
00085     =========   
00086 
00087     SIDE    (input) CHARACTER*1   
00088             = 'R': compute right eigenvectors only;   
00089             = 'L': compute left eigenvectors only;   
00090             = 'B': compute both right and left eigenvectors.   
00091 
00092     HOWMNY  (input) CHARACTER*1   
00093             = 'A': compute all right and/or left eigenvectors;   
00094             = 'B': compute all right and/or left eigenvectors, and   
00095                    backtransform them using the input matrices supplied   
00096                    in VR and/or VL;   
00097             = 'S': compute selected right and/or left eigenvectors,   
00098                    specified by the logical array SELECT.   
00099 
00100     SELECT  (input) LOGICAL array, dimension (N)   
00101             If HOWMNY='S', SELECT specifies the eigenvectors to be   
00102             computed.   
00103             If HOWMNY='A' or 'B', SELECT is not referenced.   
00104             To select the real eigenvector corresponding to the real   
00105             eigenvalue w(j), SELECT(j) must be set to .TRUE.  To select   
00106             the complex eigenvector corresponding to a complex conjugate   
00107             pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must   
00108             be set to .TRUE..   
00109 
00110     N       (input) INTEGER   
00111             The order of the matrices A and B.  N >= 0.   
00112 
00113     A       (input) DOUBLE PRECISION array, dimension (LDA,N)   
00114             The upper quasi-triangular matrix A.   
00115 
00116     LDA     (input) INTEGER   
00117             The leading dimension of array A.  LDA >= max(1, N).   
00118 
00119     B       (input) DOUBLE PRECISION array, dimension (LDB,N)   
00120             The upper triangular matrix B.  If A has a 2-by-2 diagonal   
00121             block, then the corresponding 2-by-2 block of B must be   
00122             diagonal with positive elements.   
00123 
00124     LDB     (input) INTEGER   
00125             The leading dimension of array B.  LDB >= max(1,N).   
00126 
00127     VL      (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)   
00128             On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must   
00129             contain an N-by-N matrix Q (usually the orthogonal matrix Q   
00130             of left Schur vectors returned by DHGEQZ).   
00131             On exit, if SIDE = 'L' or 'B', VL contains:   
00132             if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B);   
00133             if HOWMNY = 'B', the matrix Q*Y;   
00134             if HOWMNY = 'S', the left eigenvectors of (A,B) specified by   
00135                         SELECT, stored consecutively in the columns of   
00136                         VL, in the same order as their eigenvalues.   
00137             If SIDE = 'R', VL is not referenced.   
00138 
00139             A complex eigenvector corresponding to a complex eigenvalue   
00140             is stored in two consecutive columns, the first holding the   
00141             real part, and the second the imaginary part.   
00142 
00143     LDVL    (input) INTEGER   
00144             The leading dimension of array VL.   
00145             LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.   
00146 
00147     VR      (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)   
00148             On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must   
00149             contain an N-by-N matrix Q (usually the orthogonal matrix Z   
00150             of right Schur vectors returned by DHGEQZ).   
00151             On exit, if SIDE = 'R' or 'B', VR contains:   
00152             if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B);   
00153             if HOWMNY = 'B', the matrix Z*X;   
00154             if HOWMNY = 'S', the right eigenvectors of (A,B) specified by   
00155                         SELECT, stored consecutively in the columns of   
00156                         VR, in the same order as their eigenvalues.   
00157             If SIDE = 'L', VR is not referenced.   
00158 
00159             A complex eigenvector corresponding to a complex eigenvalue   
00160             is stored in two consecutive columns, the first holding the   
00161             real part and the second the imaginary part.   
00162 
00163     LDVR    (input) INTEGER   
00164             The leading dimension of the array VR.   
00165             LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.   
00166 
00167     MM      (input) INTEGER   
00168             The number of columns in the arrays VL and/or VR. MM >= M.   
00169 
00170     M       (output) INTEGER   
00171             The number of columns in the arrays VL and/or VR actually   
00172             used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M   
00173             is set to N.  Each selected real eigenvector occupies one   
00174             column and each selected complex eigenvector occupies two   
00175             columns.   
00176 
00177     WORK    (workspace) DOUBLE PRECISION array, dimension (6*N)   
00178 
00179     INFO    (output) INTEGER   
00180             = 0:  successful exit.   
00181             < 0:  if INFO = -i, the i-th argument had an illegal value.   
00182             > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex   
00183                   eigenvalue.   
00184 
00185     Further Details   
00186     ===============   
00187 
00188     Allocation of workspace:   
00189     ---------- -- ---------   
00190 
00191        WORK( j ) = 1-norm of j-th column of A, above the diagonal   
00192        WORK( N+j ) = 1-norm of j-th column of B, above the diagonal   
00193        WORK( 2*N+1:3*N ) = real part of eigenvector   
00194        WORK( 3*N+1:4*N ) = imaginary part of eigenvector   
00195        WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector   
00196        WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector   
00197 
00198     Rowwise vs. columnwise solution methods:   
00199     ------- --  ---------- -------- -------   
00200 
00201     Finding a generalized eigenvector consists basically of solving the   
00202     singular triangular system   
00203 
00204      (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)   
00205 
00206     Consider finding the i-th right eigenvector (assume all eigenvalues   
00207     are real). The equation to be solved is:   
00208          n                   i   
00209     0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1   
00210         k=j                 k=j   
00211 
00212     where  C = (A - w B)  (The components v(i+1:n) are 0.)   
00213 
00214     The "rowwise" method is:   
00215 
00216     (1)  v(i) := 1   
00217     for j = i-1,. . .,1:   
00218                             i   
00219         (2) compute  s = - sum C(j,k) v(k)   and   
00220                           k=j+1   
00221 
00222         (3) v(j) := s / C(j,j)   
00223 
00224     Step 2 is sometimes called the "dot product" step, since it is an   
00225     inner product between the j-th row and the portion of the eigenvector   
00226     that has been computed so far.   
00227 
00228     The "columnwise" method consists basically in doing the sums   
00229     for all the rows in parallel.  As each v(j) is computed, the   
00230     contribution of v(j) times the j-th column of C is added to the   
00231     partial sums.  Since FORTRAN arrays are stored columnwise, this has   
00232     the advantage that at each step, the elements of C that are accessed   
00233     are adjacent to one another, whereas with the rowwise method, the   
00234     elements accessed at a step are spaced LDA (and LDB) words apart.   
00235 
00236     When finding left eigenvectors, the matrix in question is the   
00237     transpose of the one in storage, so the rowwise method then   
00238     actually accesses columns of A and B at each step, and so is the   
00239     preferred method.   
00240 
00241     =====================================================================   
00242 
00243 
00244        Decode and Test the input parameters   
00245 
00246        Parameter adjustments */
00247     /* Table of constant values */
00248      logical c_true = TRUE_;
00249      integer c__2 = 2;
00250      Treal c_b35 = 1.;
00251      integer c__1 = 1;
00252      Treal c_b37 = 0.;
00253      logical c_false = FALSE_;
00254     
00255     /* System generated locals */
00256     integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1, 
00257             vr_offset, i__1, i__2, i__3, i__4, i__5;
00258     Treal d__1, d__2, d__3, d__4, d__5, d__6;
00259     /* Local variables */
00260      integer ibeg, ieig, iend;
00261      Treal dmin__, temp, suma[4]        /* was [2][2] */, sumb[4]       
00262             /* was [2][2] */, xmax;
00263      Treal cim2a, cim2b, cre2a, cre2b, temp2, bdiag[2];
00264      integer i__, j;
00265      Treal acoef, scale;
00266      logical ilall;
00267      integer iside;
00268      Treal sbeta;
00269      logical il2by2;
00270      integer iinfo;
00271      Treal small;
00272      logical compl_AAAA;
00273      Treal anorm, bnorm;
00274      logical compr;
00275      Treal temp2i;
00276      Treal temp2r;
00277      integer ja;
00278      logical ilabad, ilbbad;
00279      integer jc, je, na;
00280      Treal acoefa, bcoefa, cimaga, cimagb;
00281      logical ilback;
00282      integer im;
00283      Treal bcoefi, ascale, bscale, creala;
00284      integer jr;
00285      Treal crealb;
00286      Treal bcoefr;
00287      integer jw, nw;
00288      Treal salfar, safmin;
00289      Treal xscale, bignum;
00290      logical ilcomp, ilcplx;
00291      integer ihwmny;
00292      Treal big;
00293      logical lsa, lsb;
00294      Treal ulp, sum[4]  /* was [2][2] */;
00295 #define suma_ref(a_1,a_2) suma[(a_2)*2 + a_1 - 3]
00296 #define sumb_ref(a_1,a_2) sumb[(a_2)*2 + a_1 - 3]
00297 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00298 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00299 #define vl_ref(a_1,a_2) vl[(a_2)*vl_dim1 + a_1]
00300 #define vr_ref(a_1,a_2) vr[(a_2)*vr_dim1 + a_1]
00301 #define sum_ref(a_1,a_2) sum[(a_2)*2 + a_1 - 3]
00302 
00303 
00304     --select;
00305     a_dim1 = *lda;
00306     a_offset = 1 + a_dim1 * 1;
00307     a -= a_offset;
00308     b_dim1 = *ldb;
00309     b_offset = 1 + b_dim1 * 1;
00310     b -= b_offset;
00311     vl_dim1 = *ldvl;
00312     vl_offset = 1 + vl_dim1 * 1;
00313     vl -= vl_offset;
00314     vr_dim1 = *ldvr;
00315     vr_offset = 1 + vr_dim1 * 1;
00316     vr -= vr_offset;
00317     --work;
00318 
00319     /* Initialization added by Elias to get rid of compiler warnings. */
00320     ilback = 0;
00321     /* Function Body */
00322     if (template_blas_lsame(howmny, "A")) {
00323         ihwmny = 1;
00324         ilall = TRUE_;
00325         ilback = FALSE_;
00326     } else if (template_blas_lsame(howmny, "S")) {
00327         ihwmny = 2;
00328         ilall = FALSE_;
00329         ilback = FALSE_;
00330     } else if (template_blas_lsame(howmny, "B") || template_blas_lsame(howmny, 
00331             "T")) {
00332         ihwmny = 3;
00333         ilall = TRUE_;
00334         ilback = TRUE_;
00335     } else {
00336         ihwmny = -1;
00337         ilall = TRUE_;
00338     }
00339 
00340     if (template_blas_lsame(side, "R")) {
00341         iside = 1;
00342         compl_AAAA = FALSE_;
00343         compr = TRUE_;
00344     } else if (template_blas_lsame(side, "L")) {
00345         iside = 2;
00346         compl_AAAA = TRUE_;
00347         compr = FALSE_;
00348     } else if (template_blas_lsame(side, "B")) {
00349         iside = 3;
00350         compl_AAAA = TRUE_;
00351         compr = TRUE_;
00352     } else {
00353         iside = -1;
00354     }
00355 
00356     *info = 0;
00357     if (iside < 0) {
00358         *info = -1;
00359     } else if (ihwmny < 0) {
00360         *info = -2;
00361     } else if (*n < 0) {
00362         *info = -4;
00363     } else if (*lda < maxMACRO(1,*n)) {
00364         *info = -6;
00365     } else if (*ldb < maxMACRO(1,*n)) {
00366         *info = -8;
00367     }
00368     if (*info != 0) {
00369         i__1 = -(*info);
00370         template_blas_erbla("TGEVC ", &i__1);
00371         return 0;
00372     }
00373 
00374 /*     Count the number of eigenvectors to be computed */
00375 
00376     if (! ilall) {
00377         im = 0;
00378         ilcplx = FALSE_;
00379         i__1 = *n;
00380         for (j = 1; j <= i__1; ++j) {
00381             if (ilcplx) {
00382                 ilcplx = FALSE_;
00383                 goto L10;
00384             }
00385             if (j < *n) {
00386                 if (a_ref(j + 1, j) != 0.) {
00387                     ilcplx = TRUE_;
00388                 }
00389             }
00390             if (ilcplx) {
00391                 if (select[j] || select[j + 1]) {
00392                     im += 2;
00393                 }
00394             } else {
00395                 if (select[j]) {
00396                     ++im;
00397                 }
00398             }
00399 L10:
00400             ;
00401         }
00402     } else {
00403         im = *n;
00404     }
00405 
00406 /*     Check 2-by-2 diagonal blocks of A, B */
00407 
00408     ilabad = FALSE_;
00409     ilbbad = FALSE_;
00410     i__1 = *n - 1;
00411     for (j = 1; j <= i__1; ++j) {
00412         if (a_ref(j + 1, j) != 0.) {
00413             if (b_ref(j, j) == 0. || b_ref(j + 1, j + 1) == 0. || b_ref(j, j 
00414                     + 1) != 0.) {
00415                 ilbbad = TRUE_;
00416             }
00417             if (j < *n - 1) {
00418                 if (a_ref(j + 2, j + 1) != 0.) {
00419                     ilabad = TRUE_;
00420                 }
00421             }
00422         }
00423 /* L20: */
00424     }
00425 
00426     if (ilabad) {
00427         *info = -5;
00428     } else if (ilbbad) {
00429         *info = -7;
00430     } else if ( ( compl_AAAA && *ldvl < *n ) || *ldvl < 1) {
00431         *info = -10;
00432     } else if ( ( compr && *ldvr < *n ) || *ldvr < 1) {
00433         *info = -12;
00434     } else if (*mm < im) {
00435         *info = -13;
00436     }
00437     if (*info != 0) {
00438         i__1 = -(*info);
00439         template_blas_erbla("TGEVC ", &i__1);
00440         return 0;
00441     }
00442 
00443 /*     Quick return if possible */
00444 
00445     *m = im;
00446     if (*n == 0) {
00447         return 0;
00448     }
00449 
00450 /*     Machine Constants */
00451 
00452     safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00453     big = 1. / safmin;
00454     template_lapack_labad(&safmin, &big);
00455     ulp = template_lapack_lamch("Epsilon", (Treal)0) * template_lapack_lamch("Base", (Treal)0);
00456     small = safmin * *n / ulp;
00457     big = 1. / small;
00458     bignum = 1. / (safmin * *n);
00459 
00460 /*     Compute the 1-norm of each column of the strictly upper triangular   
00461        part (i.e., excluding all elements belonging to the diagonal   
00462        blocks) of A and B to check for possible overflow in the   
00463        triangular solver. */
00464 
00465     anorm = (d__1 = a_ref(1, 1), absMACRO(d__1));
00466     if (*n > 1) {
00467         anorm += (d__1 = a_ref(2, 1), absMACRO(d__1));
00468     }
00469     bnorm = (d__1 = b_ref(1, 1), absMACRO(d__1));
00470     work[1] = 0.;
00471     work[*n + 1] = 0.;
00472 
00473     i__1 = *n;
00474     for (j = 2; j <= i__1; ++j) {
00475         temp = 0.;
00476         temp2 = 0.;
00477         if (a_ref(j, j - 1) == 0.) {
00478             iend = j - 1;
00479         } else {
00480             iend = j - 2;
00481         }
00482         i__2 = iend;
00483         for (i__ = 1; i__ <= i__2; ++i__) {
00484             temp += (d__1 = a_ref(i__, j), absMACRO(d__1));
00485             temp2 += (d__1 = b_ref(i__, j), absMACRO(d__1));
00486 /* L30: */
00487         }
00488         work[j] = temp;
00489         work[*n + j] = temp2;
00490 /* Computing MIN */
00491         i__3 = j + 1;
00492         i__2 = minMACRO(i__3,*n);
00493         for (i__ = iend + 1; i__ <= i__2; ++i__) {
00494             temp += (d__1 = a_ref(i__, j), absMACRO(d__1));
00495             temp2 += (d__1 = b_ref(i__, j), absMACRO(d__1));
00496 /* L40: */
00497         }
00498         anorm = maxMACRO(anorm,temp);
00499         bnorm = maxMACRO(bnorm,temp2);
00500 /* L50: */
00501     }
00502 
00503     ascale = 1. / maxMACRO(anorm,safmin);
00504     bscale = 1. / maxMACRO(bnorm,safmin);
00505 
00506 /*     Left eigenvectors */
00507 
00508     if (compl_AAAA) {
00509         ieig = 0;
00510 
00511 /*        Main loop over eigenvalues */
00512 
00513         ilcplx = FALSE_;
00514         i__1 = *n;
00515         for (je = 1; je <= i__1; ++je) {
00516 
00517 /*           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or   
00518              (b) this would be the second of a complex pair.   
00519              Check for complex eigenvalue, so as to be sure of which   
00520              entry(-ies) of SELECT to look at. */
00521 
00522             if (ilcplx) {
00523                 ilcplx = FALSE_;
00524                 goto L220;
00525             }
00526             nw = 1;
00527             if (je < *n) {
00528                 if (a_ref(je + 1, je) != 0.) {
00529                     ilcplx = TRUE_;
00530                     nw = 2;
00531                 }
00532             }
00533             if (ilall) {
00534                 ilcomp = TRUE_;
00535             } else if (ilcplx) {
00536                 ilcomp = select[je] || select[je + 1];
00537             } else {
00538                 ilcomp = select[je];
00539             }
00540             if (! ilcomp) {
00541                 goto L220;
00542             }
00543 
00544 /*           Decide if (a) singular pencil, (b) real eigenvalue, or   
00545              (c) complex eigenvalue. */
00546 
00547             if (! ilcplx) {
00548                 if ((d__1 = a_ref(je, je), absMACRO(d__1)) <= safmin && (d__2 = 
00549                         b_ref(je, je), absMACRO(d__2)) <= safmin) {
00550 
00551 /*                 Singular matrix pencil -- return unit eigenvector */
00552 
00553                     ++ieig;
00554                     i__2 = *n;
00555                     for (jr = 1; jr <= i__2; ++jr) {
00556                         vl_ref(jr, ieig) = 0.;
00557 /* L60: */
00558                     }
00559                     vl_ref(ieig, ieig) = 1.;
00560                     goto L220;
00561                 }
00562             }
00563 
00564 /*           Clear vector */
00565 
00566             i__2 = nw * *n;
00567             for (jr = 1; jr <= i__2; ++jr) {
00568                 work[(*n << 1) + jr] = 0.;
00569 /* L70: */
00570             }
00571 /*                                                 T   
00572              Compute coefficients in  ( a A - b B )  y = 0   
00573                 a  is  ACOEF   
00574                 b  is  BCOEFR + i*BCOEFI */
00575 
00576             if (! ilcplx) {
00577 
00578 /*              Real eigenvalue   
00579 
00580    Computing MAX */
00581                 d__3 = (d__1 = a_ref(je, je), absMACRO(d__1)) * ascale, d__4 = (
00582                         d__2 = b_ref(je, je), absMACRO(d__2)) * bscale, d__3 = maxMACRO(
00583                         d__3,d__4);
00584                 temp = 1. / maxMACRO(d__3,safmin);
00585                 salfar = temp * a_ref(je, je) * ascale;
00586                 sbeta = temp * b_ref(je, je) * bscale;
00587                 acoef = sbeta * ascale;
00588                 bcoefr = salfar * bscale;
00589                 bcoefi = 0.;
00590 
00591 /*              Scale to avoid underflow */
00592 
00593                 scale = 1.;
00594                 lsa = absMACRO(sbeta) >= safmin && absMACRO(acoef) < small;
00595                 lsb = absMACRO(salfar) >= safmin && absMACRO(bcoefr) < small;
00596                 if (lsa) {
00597                     scale = small / absMACRO(sbeta) * minMACRO(anorm,big);
00598                 }
00599                 if (lsb) {
00600 /* Computing MAX */
00601                     d__1 = scale, d__2 = small / absMACRO(salfar) * minMACRO(bnorm,big);
00602                     scale = maxMACRO(d__1,d__2);
00603                 }
00604                 if (lsa || lsb) {
00605 /* Computing MIN   
00606    Computing MAX */
00607                     d__3 = 1., d__4 = absMACRO(acoef), d__3 = maxMACRO(d__3,d__4), d__4 
00608                             = absMACRO(bcoefr);
00609                     d__1 = scale, d__2 = 1. / (safmin * maxMACRO(d__3,d__4));
00610                     scale = minMACRO(d__1,d__2);
00611                     if (lsa) {
00612                         acoef = ascale * (scale * sbeta);
00613                     } else {
00614                         acoef = scale * acoef;
00615                     }
00616                     if (lsb) {
00617                         bcoefr = bscale * (scale * salfar);
00618                     } else {
00619                         bcoefr = scale * bcoefr;
00620                     }
00621                 }
00622                 acoefa = absMACRO(acoef);
00623                 bcoefa = absMACRO(bcoefr);
00624 
00625 /*              First component is 1 */
00626 
00627                 work[(*n << 1) + je] = 1.;
00628                 xmax = 1.;
00629             } else {
00630 
00631 /*              Complex eigenvalue */
00632 
00633                 d__1 = safmin * 100.;
00634                 template_lapack_lag2(&a_ref(je, je), lda, &b_ref(je, je), ldb, &d__1, &
00635                         acoef, &temp, &bcoefr, &temp2, &bcoefi);
00636                 bcoefi = -bcoefi;
00637                 if (bcoefi == 0.) {
00638                     *info = je;
00639                     return 0;
00640                 }
00641 
00642 /*              Scale to avoid over/underflow */
00643 
00644                 acoefa = absMACRO(acoef);
00645                 bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi);
00646                 scale = 1.;
00647                 if (acoefa * ulp < safmin && acoefa >= safmin) {
00648                     scale = safmin / ulp / acoefa;
00649                 }
00650                 if (bcoefa * ulp < safmin && bcoefa >= safmin) {
00651 /* Computing MAX */
00652                     d__1 = scale, d__2 = safmin / ulp / bcoefa;
00653                     scale = maxMACRO(d__1,d__2);
00654                 }
00655                 if (safmin * acoefa > ascale) {
00656                     scale = ascale / (safmin * acoefa);
00657                 }
00658                 if (safmin * bcoefa > bscale) {
00659 /* Computing MIN */
00660                     d__1 = scale, d__2 = bscale / (safmin * bcoefa);
00661                     scale = minMACRO(d__1,d__2);
00662                 }
00663                 if (scale != 1.) {
00664                     acoef = scale * acoef;
00665                     acoefa = absMACRO(acoef);
00666                     bcoefr = scale * bcoefr;
00667                     bcoefi = scale * bcoefi;
00668                     bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi);
00669                 }
00670 
00671 /*              Compute first two components of eigenvector */
00672 
00673                 temp = acoef * a_ref(je + 1, je);
00674                 temp2r = acoef * a_ref(je, je) - bcoefr * b_ref(je, je);
00675                 temp2i = -bcoefi * b_ref(je, je);
00676                 if (absMACRO(temp) > absMACRO(temp2r) + absMACRO(temp2i)) {
00677                     work[(*n << 1) + je] = 1.;
00678                     work[*n * 3 + je] = 0.;
00679                     work[(*n << 1) + je + 1] = -temp2r / temp;
00680                     work[*n * 3 + je + 1] = -temp2i / temp;
00681                 } else {
00682                     work[(*n << 1) + je + 1] = 1.;
00683                     work[*n * 3 + je + 1] = 0.;
00684                     temp = acoef * a_ref(je, je + 1);
00685                     work[(*n << 1) + je] = (bcoefr * b_ref(je + 1, je + 1) - 
00686                             acoef * a_ref(je + 1, je + 1)) / temp;
00687                     work[*n * 3 + je] = bcoefi * b_ref(je + 1, je + 1) / temp;
00688                 }
00689 /* Computing MAX */
00690                 d__5 = (d__1 = work[(*n << 1) + je], absMACRO(d__1)) + (d__2 = 
00691                         work[*n * 3 + je], absMACRO(d__2)), d__6 = (d__3 = work[(*
00692                         n << 1) + je + 1], absMACRO(d__3)) + (d__4 = work[*n * 3 + 
00693                         je + 1], absMACRO(d__4));
00694                 xmax = maxMACRO(d__5,d__6);
00695             }
00696 
00697 /* Computing MAX */
00698             d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 = 
00699                     maxMACRO(d__1,d__2);
00700             dmin__ = maxMACRO(d__1,safmin);
00701 
00702 /*                                           T   
00703              Triangular solve of  (a A - b B)  y = 0   
00704 
00705                                      T   
00706              (rowwise in  (a A - b B) , or columnwise in (a A - b B) ) */
00707 
00708             il2by2 = FALSE_;
00709 
00710             i__2 = *n;
00711             for (j = je + nw; j <= i__2; ++j) {
00712                 if (il2by2) {
00713                     il2by2 = FALSE_;
00714                     goto L160;
00715                 }
00716 
00717                 na = 1;
00718                 bdiag[0] = b_ref(j, j);
00719                 if (j < *n) {
00720                     if (a_ref(j + 1, j) != 0.) {
00721                         il2by2 = TRUE_;
00722                         bdiag[1] = b_ref(j + 1, j + 1);
00723                         na = 2;
00724                     }
00725                 }
00726 
00727 /*              Check whether scaling is necessary for dot products */
00728 
00729                 xscale = 1. / maxMACRO(1.,xmax);
00730 /* Computing MAX */
00731                 d__1 = work[j], d__2 = work[*n + j], d__1 = maxMACRO(d__1,d__2), 
00732                         d__2 = acoefa * work[j] + bcoefa * work[*n + j];
00733                 temp = maxMACRO(d__1,d__2);
00734                 if (il2by2) {
00735 /* Computing MAX */
00736                     d__1 = temp, d__2 = work[j + 1], d__1 = maxMACRO(d__1,d__2), 
00737                             d__2 = work[*n + j + 1], d__1 = maxMACRO(d__1,d__2), 
00738                             d__2 = acoefa * work[j + 1] + bcoefa * work[*n + 
00739                             j + 1];
00740                     temp = maxMACRO(d__1,d__2);
00741                 }
00742                 if (temp > bignum * xscale) {
00743                     i__3 = nw - 1;
00744                     for (jw = 0; jw <= i__3; ++jw) {
00745                         i__4 = j - 1;
00746                         for (jr = je; jr <= i__4; ++jr) {
00747                             work[(jw + 2) * *n + jr] = xscale * work[(jw + 2) 
00748                                     * *n + jr];
00749 /* L80: */
00750                         }
00751 /* L90: */
00752                     }
00753                     xmax *= xscale;
00754                 }
00755 
00756 /*              Compute dot products   
00757 
00758                       j-1   
00759                 SUM = sum  conjg( a*A(k,j) - b*B(k,j) )*x(k)   
00760                       k=je   
00761 
00762                 To reduce the op count, this is done as   
00763 
00764                 _        j-1                  _        j-1   
00765                 a*conjg( sum  A(k,j)*x(k) ) - b*conjg( sum  B(k,j)*x(k) )   
00766                          k=je                          k=je   
00767 
00768                 which may cause underflow problems if A or B are close   
00769                 to underflow.  (E.g., less than SMALL.)   
00770 
00771 
00772                 A series of compiler directives to defeat vectorization   
00773                 for the next loop   
00774 
00775    $PL$ CMCHAR=' '   
00776    DIR$          NEXTSCALAR   
00777    $DIR          SCALAR   
00778    DIR$          NEXT SCALAR   
00779    VD$L          NOVECTOR   
00780    DEC$          NOVECTOR   
00781    VD$           NOVECTOR   
00782    VDIR          NOVECTOR   
00783    VOCL          LOOP,SCALAR   
00784    IBM           PREFER SCALAR   
00785    $PL$ CMCHAR='*' */
00786 
00787                 i__3 = nw;
00788                 for (jw = 1; jw <= i__3; ++jw) {
00789 
00790 /* $PL$ CMCHAR=' '   
00791    DIR$             NEXTSCALAR   
00792    $DIR             SCALAR   
00793    DIR$             NEXT SCALAR   
00794    VD$L             NOVECTOR   
00795    DEC$             NOVECTOR   
00796    VD$              NOVECTOR   
00797    VDIR             NOVECTOR   
00798    VOCL             LOOP,SCALAR   
00799    IBM              PREFER SCALAR   
00800    $PL$ CMCHAR='*' */
00801 
00802                     i__4 = na;
00803                     for (ja = 1; ja <= i__4; ++ja) {
00804                         suma_ref(ja, jw) = 0.;
00805                         sumb_ref(ja, jw) = 0.;
00806 
00807                         i__5 = j - 1;
00808                         for (jr = je; jr <= i__5; ++jr) {
00809                             suma_ref(ja, jw) = suma_ref(ja, jw) + a_ref(jr, j 
00810                                     + ja - 1) * work[(jw + 1) * *n + jr];
00811                             sumb_ref(ja, jw) = sumb_ref(ja, jw) + b_ref(jr, j 
00812                                     + ja - 1) * work[(jw + 1) * *n + jr];
00813 /* L100: */
00814                         }
00815 /* L110: */
00816                     }
00817 /* L120: */
00818                 }
00819 
00820 /* $PL$ CMCHAR=' '   
00821    DIR$          NEXTSCALAR   
00822    $DIR          SCALAR   
00823    DIR$          NEXT SCALAR   
00824    VD$L          NOVECTOR   
00825    DEC$          NOVECTOR   
00826    VD$           NOVECTOR   
00827    VDIR          NOVECTOR   
00828    VOCL          LOOP,SCALAR   
00829    IBM           PREFER SCALAR   
00830    $PL$ CMCHAR='*' */
00831 
00832                 i__3 = na;
00833                 for (ja = 1; ja <= i__3; ++ja) {
00834                     if (ilcplx) {
00835                         sum_ref(ja, 1) = -acoef * suma_ref(ja, 1) + bcoefr * 
00836                                 sumb_ref(ja, 1) - bcoefi * sumb_ref(ja, 2);
00837                         sum_ref(ja, 2) = -acoef * suma_ref(ja, 2) + bcoefr * 
00838                                 sumb_ref(ja, 2) + bcoefi * sumb_ref(ja, 1);
00839                     } else {
00840                         sum_ref(ja, 1) = -acoef * suma_ref(ja, 1) + bcoefr * 
00841                                 sumb_ref(ja, 1);
00842                     }
00843 /* L130: */
00844                 }
00845 
00846 /*                                  T   
00847                 Solve  ( a A - b B )  y = SUM(,)   
00848                 with scaling and perturbation of the denominator */
00849 
00850                 template_lapack_laln2(&c_true, &na, &nw, &dmin__, &acoef, &a_ref(j, j), lda,
00851                          bdiag, &bdiag[1], sum, &c__2, &bcoefr, &bcoefi, &
00852                         work[(*n << 1) + j], n, &scale, &temp, &iinfo);
00853                 if (scale < 1.) {
00854                     i__3 = nw - 1;
00855                     for (jw = 0; jw <= i__3; ++jw) {
00856                         i__4 = j - 1;
00857                         for (jr = je; jr <= i__4; ++jr) {
00858                             work[(jw + 2) * *n + jr] = scale * work[(jw + 2) *
00859                                      *n + jr];
00860 /* L140: */
00861                         }
00862 /* L150: */
00863                     }
00864                     xmax = scale * xmax;
00865                 }
00866                 xmax = maxMACRO(xmax,temp);
00867 L160:
00868                 ;
00869             }
00870 
00871 /*           Copy eigenvector to VL, back transforming if   
00872              HOWMNY='B'. */
00873 
00874             ++ieig;
00875             if (ilback) {
00876                 i__2 = nw - 1;
00877                 for (jw = 0; jw <= i__2; ++jw) {
00878                     i__3 = *n + 1 - je;
00879                     template_blas_gemv("N", n, &i__3, &c_b35, &vl_ref(1, je), ldvl, &work[
00880                             (jw + 2) * *n + je], &c__1, &c_b37, &work[(jw + 4)
00881                              * *n + 1], &c__1);
00882 /* L170: */
00883                 }
00884                 template_lapack_lacpy(" ", n, &nw, &work[(*n << 2) + 1], n, &vl_ref(1, je), 
00885                         ldvl);
00886                 ibeg = 1;
00887             } else {
00888                 template_lapack_lacpy(" ", n, &nw, &work[(*n << 1) + 1], n, &vl_ref(1, ieig)
00889                         , ldvl);
00890                 ibeg = je;
00891             }
00892 
00893 /*           Scale eigenvector */
00894 
00895             xmax = 0.;
00896             if (ilcplx) {
00897                 i__2 = *n;
00898                 for (j = ibeg; j <= i__2; ++j) {
00899 /* Computing MAX */
00900                     d__3 = xmax, d__4 = (d__1 = vl_ref(j, ieig), absMACRO(d__1)) + 
00901                             (d__2 = vl_ref(j, ieig + 1), absMACRO(d__2));
00902                     xmax = maxMACRO(d__3,d__4);
00903 /* L180: */
00904                 }
00905             } else {
00906                 i__2 = *n;
00907                 for (j = ibeg; j <= i__2; ++j) {
00908 /* Computing MAX */
00909                     d__2 = xmax, d__3 = (d__1 = vl_ref(j, ieig), absMACRO(d__1));
00910                     xmax = maxMACRO(d__2,d__3);
00911 /* L190: */
00912                 }
00913             }
00914 
00915             if (xmax > safmin) {
00916                 xscale = 1. / xmax;
00917 
00918                 i__2 = nw - 1;
00919                 for (jw = 0; jw <= i__2; ++jw) {
00920                     i__3 = *n;
00921                     for (jr = ibeg; jr <= i__3; ++jr) {
00922                         vl_ref(jr, ieig + jw) = xscale * vl_ref(jr, ieig + jw)
00923                                 ;
00924 /* L200: */
00925                     }
00926 /* L210: */
00927                 }
00928             }
00929             ieig = ieig + nw - 1;
00930 
00931 L220:
00932             ;
00933         }
00934     }
00935 
00936 /*     Right eigenvectors */
00937 
00938     if (compr) {
00939         ieig = im + 1;
00940 
00941 /*        Main loop over eigenvalues */
00942 
00943         ilcplx = FALSE_;
00944         for (je = *n; je >= 1; --je) {
00945 
00946 /*           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or   
00947              (b) this would be the second of a complex pair.   
00948              Check for complex eigenvalue, so as to be sure of which   
00949              entry(-ies) of SELECT to look at -- if complex, SELECT(JE)   
00950              or SELECT(JE-1).   
00951              If this is a complex pair, the 2-by-2 diagonal block   
00952              corresponding to the eigenvalue is in rows/columns JE-1:JE */
00953 
00954             if (ilcplx) {
00955                 ilcplx = FALSE_;
00956                 goto L500;
00957             }
00958             nw = 1;
00959             if (je > 1) {
00960                 if (a_ref(je, je - 1) != 0.) {
00961                     ilcplx = TRUE_;
00962                     nw = 2;
00963                 }
00964             }
00965             if (ilall) {
00966                 ilcomp = TRUE_;
00967             } else if (ilcplx) {
00968                 ilcomp = select[je] || select[je - 1];
00969             } else {
00970                 ilcomp = select[je];
00971             }
00972             if (! ilcomp) {
00973                 goto L500;
00974             }
00975 
00976 /*           Decide if (a) singular pencil, (b) real eigenvalue, or   
00977              (c) complex eigenvalue. */
00978 
00979             if (! ilcplx) {
00980                 if ((d__1 = a_ref(je, je), absMACRO(d__1)) <= safmin && (d__2 = 
00981                         b_ref(je, je), absMACRO(d__2)) <= safmin) {
00982 
00983 /*                 Singular matrix pencil -- unit eigenvector */
00984 
00985                     --ieig;
00986                     i__1 = *n;
00987                     for (jr = 1; jr <= i__1; ++jr) {
00988                         vr_ref(jr, ieig) = 0.;
00989 /* L230: */
00990                     }
00991                     vr_ref(ieig, ieig) = 1.;
00992                     goto L500;
00993                 }
00994             }
00995 
00996 /*           Clear vector */
00997 
00998             i__1 = nw - 1;
00999             for (jw = 0; jw <= i__1; ++jw) {
01000                 i__2 = *n;
01001                 for (jr = 1; jr <= i__2; ++jr) {
01002                     work[(jw + 2) * *n + jr] = 0.;
01003 /* L240: */
01004                 }
01005 /* L250: */
01006             }
01007 
01008 /*           Compute coefficients in  ( a A - b B ) x = 0   
01009                 a  is  ACOEF   
01010                 b  is  BCOEFR + i*BCOEFI */
01011 
01012             if (! ilcplx) {
01013 
01014 /*              Real eigenvalue   
01015 
01016    Computing MAX */
01017                 d__3 = (d__1 = a_ref(je, je), absMACRO(d__1)) * ascale, d__4 = (
01018                         d__2 = b_ref(je, je), absMACRO(d__2)) * bscale, d__3 = maxMACRO(
01019                         d__3,d__4);
01020                 temp = 1. / maxMACRO(d__3,safmin);
01021                 salfar = temp * a_ref(je, je) * ascale;
01022                 sbeta = temp * b_ref(je, je) * bscale;
01023                 acoef = sbeta * ascale;
01024                 bcoefr = salfar * bscale;
01025                 bcoefi = 0.;
01026 
01027 /*              Scale to avoid underflow */
01028 
01029                 scale = 1.;
01030                 lsa = absMACRO(sbeta) >= safmin && absMACRO(acoef) < small;
01031                 lsb = absMACRO(salfar) >= safmin && absMACRO(bcoefr) < small;
01032                 if (lsa) {
01033                     scale = small / absMACRO(sbeta) * minMACRO(anorm,big);
01034                 }
01035                 if (lsb) {
01036 /* Computing MAX */
01037                     d__1 = scale, d__2 = small / absMACRO(salfar) * minMACRO(bnorm,big);
01038                     scale = maxMACRO(d__1,d__2);
01039                 }
01040                 if (lsa || lsb) {
01041 /* Computing MIN   
01042    Computing MAX */
01043                     d__3 = 1., d__4 = absMACRO(acoef), d__3 = maxMACRO(d__3,d__4), d__4 
01044                             = absMACRO(bcoefr);
01045                     d__1 = scale, d__2 = 1. / (safmin * maxMACRO(d__3,d__4));
01046                     scale = minMACRO(d__1,d__2);
01047                     if (lsa) {
01048                         acoef = ascale * (scale * sbeta);
01049                     } else {
01050                         acoef = scale * acoef;
01051                     }
01052                     if (lsb) {
01053                         bcoefr = bscale * (scale * salfar);
01054                     } else {
01055                         bcoefr = scale * bcoefr;
01056                     }
01057                 }
01058                 acoefa = absMACRO(acoef);
01059                 bcoefa = absMACRO(bcoefr);
01060 
01061 /*              First component is 1 */
01062 
01063                 work[(*n << 1) + je] = 1.;
01064                 xmax = 1.;
01065 
01066 /*              Compute contribution from column JE of A and B to sum   
01067                 (See "Further Details", above.) */
01068 
01069                 i__1 = je - 1;
01070                 for (jr = 1; jr <= i__1; ++jr) {
01071                     work[(*n << 1) + jr] = bcoefr * b_ref(jr, je) - acoef * 
01072                             a_ref(jr, je);
01073 /* L260: */
01074                 }
01075             } else {
01076 
01077 /*              Complex eigenvalue */
01078 
01079                 d__1 = safmin * 100.;
01080                 template_lapack_lag2(&a_ref(je - 1, je - 1), lda, &b_ref(je - 1, je - 1), 
01081                         ldb, &d__1, &acoef, &temp, &bcoefr, &temp2, &bcoefi);
01082                 if (bcoefi == 0.) {
01083                     *info = je - 1;
01084                     return 0;
01085                 }
01086 
01087 /*              Scale to avoid over/underflow */
01088 
01089                 acoefa = absMACRO(acoef);
01090                 bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi);
01091                 scale = 1.;
01092                 if (acoefa * ulp < safmin && acoefa >= safmin) {
01093                     scale = safmin / ulp / acoefa;
01094                 }
01095                 if (bcoefa * ulp < safmin && bcoefa >= safmin) {
01096 /* Computing MAX */
01097                     d__1 = scale, d__2 = safmin / ulp / bcoefa;
01098                     scale = maxMACRO(d__1,d__2);
01099                 }
01100                 if (safmin * acoefa > ascale) {
01101                     scale = ascale / (safmin * acoefa);
01102                 }
01103                 if (safmin * bcoefa > bscale) {
01104 /* Computing MIN */
01105                     d__1 = scale, d__2 = bscale / (safmin * bcoefa);
01106                     scale = minMACRO(d__1,d__2);
01107                 }
01108                 if (scale != 1.) {
01109                     acoef = scale * acoef;
01110                     acoefa = absMACRO(acoef);
01111                     bcoefr = scale * bcoefr;
01112                     bcoefi = scale * bcoefi;
01113                     bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi);
01114                 }
01115 
01116 /*              Compute first two components of eigenvector   
01117                 and contribution to sums */
01118 
01119                 temp = acoef * a_ref(je, je - 1);
01120                 temp2r = acoef * a_ref(je, je) - bcoefr * b_ref(je, je);
01121                 temp2i = -bcoefi * b_ref(je, je);
01122                 if (absMACRO(temp) >= absMACRO(temp2r) + absMACRO(temp2i)) {
01123                     work[(*n << 1) + je] = 1.;
01124                     work[*n * 3 + je] = 0.;
01125                     work[(*n << 1) + je - 1] = -temp2r / temp;
01126                     work[*n * 3 + je - 1] = -temp2i / temp;
01127                 } else {
01128                     work[(*n << 1) + je - 1] = 1.;
01129                     work[*n * 3 + je - 1] = 0.;
01130                     temp = acoef * a_ref(je - 1, je);
01131                     work[(*n << 1) + je] = (bcoefr * b_ref(je - 1, je - 1) - 
01132                             acoef * a_ref(je - 1, je - 1)) / temp;
01133                     work[*n * 3 + je] = bcoefi * b_ref(je - 1, je - 1) / temp;
01134                 }
01135 
01136 /* Computing MAX */
01137                 d__5 = (d__1 = work[(*n << 1) + je], absMACRO(d__1)) + (d__2 = 
01138                         work[*n * 3 + je], absMACRO(d__2)), d__6 = (d__3 = work[(*
01139                         n << 1) + je - 1], absMACRO(d__3)) + (d__4 = work[*n * 3 + 
01140                         je - 1], absMACRO(d__4));
01141                 xmax = maxMACRO(d__5,d__6);
01142 
01143 /*              Compute contribution from columns JE and JE-1   
01144                 of A and B to the sums. */
01145 
01146                 creala = acoef * work[(*n << 1) + je - 1];
01147                 cimaga = acoef * work[*n * 3 + je - 1];
01148                 crealb = bcoefr * work[(*n << 1) + je - 1] - bcoefi * work[*n 
01149                         * 3 + je - 1];
01150                 cimagb = bcoefi * work[(*n << 1) + je - 1] + bcoefr * work[*n 
01151                         * 3 + je - 1];
01152                 cre2a = acoef * work[(*n << 1) + je];
01153                 cim2a = acoef * work[*n * 3 + je];
01154                 cre2b = bcoefr * work[(*n << 1) + je] - bcoefi * work[*n * 3 
01155                         + je];
01156                 cim2b = bcoefi * work[(*n << 1) + je] + bcoefr * work[*n * 3 
01157                         + je];
01158                 i__1 = je - 2;
01159                 for (jr = 1; jr <= i__1; ++jr) {
01160                     work[(*n << 1) + jr] = -creala * a_ref(jr, je - 1) + 
01161                             crealb * b_ref(jr, je - 1) - cre2a * a_ref(jr, je)
01162                              + cre2b * b_ref(jr, je);
01163                     work[*n * 3 + jr] = -cimaga * a_ref(jr, je - 1) + cimagb *
01164                              b_ref(jr, je - 1) - cim2a * a_ref(jr, je) + 
01165                             cim2b * b_ref(jr, je);
01166 /* L270: */
01167                 }
01168             }
01169 
01170 /* Computing MAX */
01171             d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 = 
01172                     maxMACRO(d__1,d__2);
01173             dmin__ = maxMACRO(d__1,safmin);
01174 
01175 /*           Columnwise triangular solve of  (a A - b B)  x = 0 */
01176 
01177             il2by2 = FALSE_;
01178             for (j = je - nw; j >= 1; --j) {
01179 
01180 /*              If a 2-by-2 block, is in position j-1:j, wait until   
01181                 next iteration to process it (when it will be j:j+1) */
01182 
01183                 if (! il2by2 && j > 1) {
01184                     if (a_ref(j, j - 1) != 0.) {
01185                         il2by2 = TRUE_;
01186                         goto L370;
01187                     }
01188                 }
01189                 bdiag[0] = b_ref(j, j);
01190                 if (il2by2) {
01191                     na = 2;
01192                     bdiag[1] = b_ref(j + 1, j + 1);
01193                 } else {
01194                     na = 1;
01195                 }
01196 
01197 /*              Compute x(j) (and x(j+1), if 2-by-2 block) */
01198 
01199                 template_lapack_laln2(&c_false, &na, &nw, &dmin__, &acoef, &a_ref(j, j), 
01200                         lda, bdiag, &bdiag[1], &work[(*n << 1) + j], n, &
01201                         bcoefr, &bcoefi, sum, &c__2, &scale, &temp, &iinfo);
01202                 if (scale < 1.) {
01203 
01204                     i__1 = nw - 1;
01205                     for (jw = 0; jw <= i__1; ++jw) {
01206                         i__2 = je;
01207                         for (jr = 1; jr <= i__2; ++jr) {
01208                             work[(jw + 2) * *n + jr] = scale * work[(jw + 2) *
01209                                      *n + jr];
01210 /* L280: */
01211                         }
01212 /* L290: */
01213                     }
01214                 }
01215 /* Computing MAX */
01216                 d__1 = scale * xmax;
01217                 xmax = maxMACRO(d__1,temp);
01218 
01219                 i__1 = nw;
01220                 for (jw = 1; jw <= i__1; ++jw) {
01221                     i__2 = na;
01222                     for (ja = 1; ja <= i__2; ++ja) {
01223                         work[(jw + 1) * *n + j + ja - 1] = sum_ref(ja, jw);
01224 /* L300: */
01225                     }
01226 /* L310: */
01227                 }
01228 
01229 /*              w = w + x(j)*(a A(*,j) - b B(*,j) ) with scaling */
01230 
01231                 if (j > 1) {
01232 
01233 /*                 Check whether scaling is necessary for sum. */
01234 
01235                     xscale = 1. / maxMACRO(1.,xmax);
01236                     temp = acoefa * work[j] + bcoefa * work[*n + j];
01237                     if (il2by2) {
01238 /* Computing MAX */
01239                         d__1 = temp, d__2 = acoefa * work[j + 1] + bcoefa * 
01240                                 work[*n + j + 1];
01241                         temp = maxMACRO(d__1,d__2);
01242                     }
01243 /* Computing MAX */
01244                     d__1 = maxMACRO(temp,acoefa);
01245                     temp = maxMACRO(d__1,bcoefa);
01246                     if (temp > bignum * xscale) {
01247 
01248                         i__1 = nw - 1;
01249                         for (jw = 0; jw <= i__1; ++jw) {
01250                             i__2 = je;
01251                             for (jr = 1; jr <= i__2; ++jr) {
01252                                 work[(jw + 2) * *n + jr] = xscale * work[(jw 
01253                                         + 2) * *n + jr];
01254 /* L320: */
01255                             }
01256 /* L330: */
01257                         }
01258                         xmax *= xscale;
01259                     }
01260 
01261 /*                 Compute the contributions of the off-diagonals of   
01262                    column j (and j+1, if 2-by-2 block) of A and B to the   
01263                    sums. */
01264 
01265 
01266                     i__1 = na;
01267                     for (ja = 1; ja <= i__1; ++ja) {
01268                         if (ilcplx) {
01269                             creala = acoef * work[(*n << 1) + j + ja - 1];
01270                             cimaga = acoef * work[*n * 3 + j + ja - 1];
01271                             crealb = bcoefr * work[(*n << 1) + j + ja - 1] - 
01272                                     bcoefi * work[*n * 3 + j + ja - 1];
01273                             cimagb = bcoefi * work[(*n << 1) + j + ja - 1] + 
01274                                     bcoefr * work[*n * 3 + j + ja - 1];
01275                             i__2 = j - 1;
01276                             for (jr = 1; jr <= i__2; ++jr) {
01277                                 work[(*n << 1) + jr] = work[(*n << 1) + jr] - 
01278                                         creala * a_ref(jr, j + ja - 1) + 
01279                                         crealb * b_ref(jr, j + ja - 1);
01280                                 work[*n * 3 + jr] = work[*n * 3 + jr] - 
01281                                         cimaga * a_ref(jr, j + ja - 1) + 
01282                                         cimagb * b_ref(jr, j + ja - 1);
01283 /* L340: */
01284                             }
01285                         } else {
01286                             creala = acoef * work[(*n << 1) + j + ja - 1];
01287                             crealb = bcoefr * work[(*n << 1) + j + ja - 1];
01288                             i__2 = j - 1;
01289                             for (jr = 1; jr <= i__2; ++jr) {
01290                                 work[(*n << 1) + jr] = work[(*n << 1) + jr] - 
01291                                         creala * a_ref(jr, j + ja - 1) + 
01292                                         crealb * b_ref(jr, j + ja - 1);
01293 /* L350: */
01294                             }
01295                         }
01296 /* L360: */
01297                     }
01298                 }
01299 
01300                 il2by2 = FALSE_;
01301 L370:
01302                 ;
01303             }
01304 
01305 /*           Copy eigenvector to VR, back transforming if   
01306              HOWMNY='B'. */
01307 
01308             ieig -= nw;
01309             if (ilback) {
01310 
01311                 i__1 = nw - 1;
01312                 for (jw = 0; jw <= i__1; ++jw) {
01313                     i__2 = *n;
01314                     for (jr = 1; jr <= i__2; ++jr) {
01315                         work[(jw + 4) * *n + jr] = work[(jw + 2) * *n + 1] * 
01316                                 vr_ref(jr, 1);
01317 /* L380: */
01318                     }
01319 
01320 /*                 A series of compiler directives to defeat   
01321                    vectorization for the next loop */
01322 
01323 
01324                     i__2 = je;
01325                     for (jc = 2; jc <= i__2; ++jc) {
01326                         i__3 = *n;
01327                         for (jr = 1; jr <= i__3; ++jr) {
01328                             work[(jw + 4) * *n + jr] += work[(jw + 2) * *n + 
01329                                     jc] * vr_ref(jr, jc);
01330 /* L390: */
01331                         }
01332 /* L400: */
01333                     }
01334 /* L410: */
01335                 }
01336 
01337                 i__1 = nw - 1;
01338                 for (jw = 0; jw <= i__1; ++jw) {
01339                     i__2 = *n;
01340                     for (jr = 1; jr <= i__2; ++jr) {
01341                         vr_ref(jr, ieig + jw) = work[(jw + 4) * *n + jr];
01342 /* L420: */
01343                     }
01344 /* L430: */
01345                 }
01346 
01347                 iend = *n;
01348             } else {
01349                 i__1 = nw - 1;
01350                 for (jw = 0; jw <= i__1; ++jw) {
01351                     i__2 = *n;
01352                     for (jr = 1; jr <= i__2; ++jr) {
01353                         vr_ref(jr, ieig + jw) = work[(jw + 2) * *n + jr];
01354 /* L440: */
01355                     }
01356 /* L450: */
01357                 }
01358 
01359                 iend = je;
01360             }
01361 
01362 /*           Scale eigenvector */
01363 
01364             xmax = 0.;
01365             if (ilcplx) {
01366                 i__1 = iend;
01367                 for (j = 1; j <= i__1; ++j) {
01368 /* Computing MAX */
01369                     d__3 = xmax, d__4 = (d__1 = vr_ref(j, ieig), absMACRO(d__1)) + 
01370                             (d__2 = vr_ref(j, ieig + 1), absMACRO(d__2));
01371                     xmax = maxMACRO(d__3,d__4);
01372 /* L460: */
01373                 }
01374             } else {
01375                 i__1 = iend;
01376                 for (j = 1; j <= i__1; ++j) {
01377 /* Computing MAX */
01378                     d__2 = xmax, d__3 = (d__1 = vr_ref(j, ieig), absMACRO(d__1));
01379                     xmax = maxMACRO(d__2,d__3);
01380 /* L470: */
01381                 }
01382             }
01383 
01384             if (xmax > safmin) {
01385                 xscale = 1. / xmax;
01386                 i__1 = nw - 1;
01387                 for (jw = 0; jw <= i__1; ++jw) {
01388                     i__2 = iend;
01389                     for (jr = 1; jr <= i__2; ++jr) {
01390                         vr_ref(jr, ieig + jw) = xscale * vr_ref(jr, ieig + jw)
01391                                 ;
01392 /* L480: */
01393                     }
01394 /* L490: */
01395                 }
01396             }
01397 L500:
01398             ;
01399         }
01400     }
01401 
01402     return 0;
01403 
01404 /*     End of DTGEVC */
01405 
01406 } /* dtgevc_ */
01407 
01408 #undef sum_ref
01409 #undef vr_ref
01410 #undef vl_ref
01411 #undef b_ref
01412 #undef a_ref
01413 #undef sumb_ref
01414 #undef suma_ref
01415 
01416 
01417 #endif

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1