template_lapack_ggev.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_GGEV_HEADER
00036 #define TEMPLATE_LAPACK_GGEV_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_ggev(const char *jobvl, const char *jobvr, const integer *n, Treal *
00041         a, const integer *lda, Treal *b, const integer *ldb, Treal *alphar, 
00042         Treal *alphai, Treal *beta, Treal *vl, const integer *ldvl, 
00043         Treal *vr, const integer *ldvr, Treal *work, const integer *lwork, 
00044         integer *info)
00045 {
00046 /*  -- LAPACK driver 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     DGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B)   
00056     the generalized eigenvalues, and optionally, the left and/or right   
00057     generalized eigenvectors.   
00058 
00059     A generalized eigenvalue for a pair of matrices (A,B) is a scalar   
00060     lambda or a ratio alpha/beta = lambda, such that A - lambda*B is   
00061     singular. It is usually represented as the pair (alpha,beta), as   
00062     there is a reasonable interpretation for beta=0, and even for both   
00063     being zero.   
00064 
00065     The right eigenvector v(j) corresponding to the eigenvalue lambda(j)   
00066     of (A,B) satisfies   
00067 
00068                      A * v(j) = lambda(j) * B * v(j).   
00069 
00070     The left eigenvector u(j) corresponding to the eigenvalue lambda(j)   
00071     of (A,B) satisfies   
00072 
00073                      u(j)**H * A  = lambda(j) * u(j)**H * B .   
00074 
00075     where u(j)**H is the conjugate-transpose of u(j).   
00076 
00077 
00078     Arguments   
00079     =========   
00080 
00081     JOBVL   (input) CHARACTER*1   
00082             = 'N':  do not compute the left generalized eigenvectors;   
00083             = 'V':  compute the left generalized eigenvectors.   
00084 
00085     JOBVR   (input) CHARACTER*1   
00086             = 'N':  do not compute the right generalized eigenvectors;   
00087             = 'V':  compute the right generalized eigenvectors.   
00088 
00089     N       (input) INTEGER   
00090             The order of the matrices A, B, VL, and VR.  N >= 0.   
00091 
00092     A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)   
00093             On entry, the matrix A in the pair (A,B).   
00094             On exit, A has been overwritten.   
00095 
00096     LDA     (input) INTEGER   
00097             The leading dimension of A.  LDA >= max(1,N).   
00098 
00099     B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)   
00100             On entry, the matrix B in the pair (A,B).   
00101             On exit, B has been overwritten.   
00102 
00103     LDB     (input) INTEGER   
00104             The leading dimension of B.  LDB >= max(1,N).   
00105 
00106     ALPHAR  (output) DOUBLE PRECISION array, dimension (N)   
00107     ALPHAI  (output) DOUBLE PRECISION array, dimension (N)   
00108     BETA    (output) DOUBLE PRECISION array, dimension (N)   
00109             On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will   
00110             be the generalized eigenvalues.  If ALPHAI(j) is zero, then   
00111             the j-th eigenvalue is real; if positive, then the j-th and   
00112             (j+1)-st eigenvalues are a complex conjugate pair, with   
00113             ALPHAI(j+1) negative.   
00114 
00115             Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)   
00116             may easily over- or underflow, and BETA(j) may even be zero.   
00117             Thus, the user should avoid naively computing the ratio   
00118             alpha/beta.  However, ALPHAR and ALPHAI will be always less   
00119             than and usually comparable with norm(A) in magnitude, and   
00120             BETA always less than and usually comparable with norm(B).   
00121 
00122     VL      (output) DOUBLE PRECISION array, dimension (LDVL,N)   
00123             If JOBVL = 'V', the left eigenvectors u(j) are stored one   
00124             after another in the columns of VL, in the same order as   
00125             their eigenvalues. If the j-th eigenvalue is real, then   
00126             u(j) = VL(:,j), the j-th column of VL. If the j-th and   
00127             (j+1)-th eigenvalues form a complex conjugate pair, then   
00128             u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).   
00129             Each eigenvector will be scaled so the largest component have   
00130             abs(real part)+abs(imag. part)=1.   
00131             Not referenced if JOBVL = 'N'.   
00132 
00133     LDVL    (input) INTEGER   
00134             The leading dimension of the matrix VL. LDVL >= 1, and   
00135             if JOBVL = 'V', LDVL >= N.   
00136 
00137     VR      (output) DOUBLE PRECISION array, dimension (LDVR,N)   
00138             If JOBVR = 'V', the right eigenvectors v(j) are stored one   
00139             after another in the columns of VR, in the same order as   
00140             their eigenvalues. If the j-th eigenvalue is real, then   
00141             v(j) = VR(:,j), the j-th column of VR. If the j-th and   
00142             (j+1)-th eigenvalues form a complex conjugate pair, then   
00143             v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).   
00144             Each eigenvector will be scaled so the largest component have   
00145             abs(real part)+abs(imag. part)=1.   
00146             Not referenced if JOBVR = 'N'.   
00147 
00148     LDVR    (input) INTEGER   
00149             The leading dimension of the matrix VR. LDVR >= 1, and   
00150             if JOBVR = 'V', LDVR >= N.   
00151 
00152     WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)   
00153             On exit, if INFO = 0, WORK(1) returns the optimal LWORK.   
00154 
00155     LWORK   (input) INTEGER   
00156             The dimension of the array WORK.  LWORK >= max(1,8*N).   
00157             For good performance, LWORK must generally be larger.   
00158 
00159             If LWORK = -1, then a workspace query is assumed; the routine   
00160             only calculates the optimal size of the WORK array, returns   
00161             this value as the first entry of the WORK array, and no error   
00162             message related to LWORK is issued by XERBLA.   
00163 
00164     INFO    (output) INTEGER   
00165             = 0:  successful exit   
00166             < 0:  if INFO = -i, the i-th argument had an illegal value.   
00167             = 1,...,N:   
00168                   The QZ iteration failed.  No eigenvectors have been   
00169                   calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)   
00170                   should be correct for j=INFO+1,...,N.   
00171             > N:  =N+1: other than QZ iteration failed in DHGEQZ.   
00172                   =N+2: error return from DTGEVC.   
00173 
00174     =====================================================================   
00175 
00176 
00177        Decode the input arguments   
00178 
00179        Parameter adjustments */
00180     /* Table of constant values */
00181      integer c__1 = 1;
00182      integer c__0 = 0;
00183      Treal c_b26 = 0.;
00184      Treal c_b27 = 1.;
00185     
00186     /* System generated locals */
00187     integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1, 
00188             vr_offset, i__1, i__2;
00189     Treal d__1, d__2, d__3, d__4;
00190     /* Local variables */
00191      Treal anrm, bnrm;
00192      integer ierr, itau;
00193      Treal temp;
00194      logical ilvl, ilvr;
00195      integer iwrk;
00196      integer ileft, icols, irows;
00197      integer jc;
00198      integer in;
00199      integer jr;
00200      logical ilascl, ilbscl;
00201      logical ldumma[1];
00202      char chtemp[1];
00203      Treal bignum;
00204      integer ijobvl, iright, ijobvr;
00205      Treal anrmto, bnrmto;
00206      integer minwrk, maxwrk;
00207      Treal smlnum;
00208      logical lquery;
00209      integer ihi, ilo;
00210      Treal eps;
00211      logical ilv;
00212 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00213 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00214 #define vl_ref(a_1,a_2) vl[(a_2)*vl_dim1 + a_1]
00215 #define vr_ref(a_1,a_2) vr[(a_2)*vr_dim1 + a_1]
00216 
00217 
00218     a_dim1 = *lda;
00219     a_offset = 1 + a_dim1 * 1;
00220     a -= a_offset;
00221     b_dim1 = *ldb;
00222     b_offset = 1 + b_dim1 * 1;
00223     b -= b_offset;
00224     --alphar;
00225     --alphai;
00226     --beta;
00227     vl_dim1 = *ldvl;
00228     vl_offset = 1 + vl_dim1 * 1;
00229     vl -= vl_offset;
00230     vr_dim1 = *ldvr;
00231     vr_offset = 1 + vr_dim1 * 1;
00232     vr -= vr_offset;
00233     --work;
00234 
00235     /* Initialization added by Elias to get rid of compiler warnings. */
00236     maxwrk = 0;
00237     /* Function Body */
00238     if (template_blas_lsame(jobvl, "N")) {
00239         ijobvl = 1;
00240         ilvl = FALSE_;
00241     } else if (template_blas_lsame(jobvl, "V")) {
00242         ijobvl = 2;
00243         ilvl = TRUE_;
00244     } else {
00245         ijobvl = -1;
00246         ilvl = FALSE_;
00247     }
00248 
00249     if (template_blas_lsame(jobvr, "N")) {
00250         ijobvr = 1;
00251         ilvr = FALSE_;
00252     } else if (template_blas_lsame(jobvr, "V")) {
00253         ijobvr = 2;
00254         ilvr = TRUE_;
00255     } else {
00256         ijobvr = -1;
00257         ilvr = FALSE_;
00258     }
00259     ilv = ilvl || ilvr;
00260 
00261 /*     Test the input arguments */
00262 
00263     *info = 0;
00264     lquery = *lwork == -1;
00265     if (ijobvl <= 0) {
00266         *info = -1;
00267     } else if (ijobvr <= 0) {
00268         *info = -2;
00269     } else if (*n < 0) {
00270         *info = -3;
00271     } else if (*lda < maxMACRO(1,*n)) {
00272         *info = -5;
00273     } else if (*ldb < maxMACRO(1,*n)) {
00274         *info = -7;
00275     } else if (*ldvl < 1 || ( ilvl && *ldvl < *n ) ) {
00276         *info = -12;
00277     } else if (*ldvr < 1 || ( ilvr && *ldvr < *n ) ) {
00278         *info = -14;
00279     }
00280 
00281 /*     Compute workspace   
00282         (Note: Comments in the code beginning "Workspace:" describe the   
00283          minimal amount of workspace needed at that point in the code,   
00284          as well as the preferred amount for good performance.   
00285          NB refers to the optimal block size for the immediately   
00286          following subroutine, as returned by ILAENV. The workspace is   
00287          computed assuming ILO = 1 and IHI = N, the worst case.) */
00288 
00289     minwrk = 1;
00290     if (*info == 0 && (*lwork >= 1 || lquery)) {
00291       maxwrk = *n * 7 + *n * template_lapack_ilaenv(&c__1, "DGEQRF", " ", n, &c__1, n, &
00292                                                     c__0, (ftnlen)6, (ftnlen)1);
00293       /* Computing MAX */
00294       i__1 = 1, i__2 = *n << 3;
00295       minwrk = maxMACRO(i__1,i__2);
00296       work[1] = (Treal) maxwrk;
00297     }
00298 
00299     if (*lwork < minwrk && ! lquery) {
00300         *info = -16;
00301     }
00302 
00303     if (*info != 0) {
00304         i__1 = -(*info);
00305         template_blas_erbla("GGEV  ", &i__1);
00306         return 0;
00307     } else if (lquery) {
00308         return 0;
00309     }
00310 
00311 /*     Quick return if possible */
00312 
00313     if (*n == 0) {
00314         return 0;
00315     }
00316 
00317 /*     Get machine constants */
00318 
00319     eps = template_lapack_lamch("P", (Treal)0);
00320     smlnum = template_lapack_lamch("S", (Treal)0);
00321     bignum = 1. / smlnum;
00322     template_lapack_labad(&smlnum, &bignum);
00323     smlnum = template_blas_sqrt(smlnum) / eps;
00324     bignum = 1. / smlnum;
00325 
00326 /*     Scale A if max element outside range [SMLNUM,BIGNUM] */
00327 
00328     anrm = template_lapack_lange("M", n, n, &a[a_offset], lda, &work[1]);
00329     ilascl = FALSE_;
00330     if (anrm > 0. && anrm < smlnum) {
00331         anrmto = smlnum;
00332         ilascl = TRUE_;
00333     } else if (anrm > bignum) {
00334         anrmto = bignum;
00335         ilascl = TRUE_;
00336     }
00337     if (ilascl) {
00338         template_lapack_lascl("G", &c__0, &c__0, &anrm, &anrmto, n, n, &a[a_offset], lda, &
00339                 ierr);
00340     }
00341 
00342 /*     Scale B if max element outside range [SMLNUM,BIGNUM] */
00343 
00344     bnrm = template_lapack_lange("M", n, n, &b[b_offset], ldb, &work[1]);
00345     ilbscl = FALSE_;
00346     if (bnrm > 0. && bnrm < smlnum) {
00347         bnrmto = smlnum;
00348         ilbscl = TRUE_;
00349     } else if (bnrm > bignum) {
00350         bnrmto = bignum;
00351         ilbscl = TRUE_;
00352     }
00353     if (ilbscl) {
00354         template_lapack_lascl("G", &c__0, &c__0, &bnrm, &bnrmto, n, n, &b[b_offset], ldb, &
00355                 ierr);
00356     }
00357 
00358 /*     Permute the matrices A, B to isolate eigenvalues if possible   
00359        (Workspace: need 6*N) */
00360 
00361     ileft = 1;
00362     iright = *n + 1;
00363     iwrk = iright + *n;
00364     template_lapack_ggbal("P", n, &a[a_offset], lda, &b[b_offset], ldb, &ilo, &ihi, &work[
00365                             ileft], &work[iright], &work[iwrk], &ierr);
00366 
00367 /*     Reduce B to triangular form (QR decomposition of B)   
00368        (Workspace: need N, prefer N*NB) */
00369 
00370     irows = ihi + 1 - ilo;
00371     if (ilv) {
00372         icols = *n + 1 - ilo;
00373     } else {
00374         icols = irows;
00375     }
00376     itau = iwrk;
00377     iwrk = itau + irows;
00378     i__1 = *lwork + 1 - iwrk;
00379     template_lapack_geqrf(&irows, &icols, &b_ref(ilo, ilo), ldb, &work[itau], &work[iwrk], &
00380             i__1, &ierr);
00381 
00382 /*     Apply the orthogonal transformation to matrix A   
00383        (Workspace: need N, prefer N*NB) */
00384 
00385     i__1 = *lwork + 1 - iwrk;
00386     /* Local char arrays added by Elias to get rid of compiler warnings. */
00387     char str_L[] = {'L', 0};
00388     char str_T[] = {'T', 0};
00389     template_lapack_ormqr(str_L, str_T, &irows, &icols, &irows, &b_ref(ilo, ilo), ldb, &work[
00390             itau], &a_ref(ilo, ilo), lda, &work[iwrk], &i__1, &ierr);
00391 
00392 /*     Initialize VL   
00393        (Workspace: need N, prefer N*NB) */
00394 
00395     if (ilvl) {
00396         template_lapack_laset("Full", n, n, &c_b26, &c_b27, &vl[vl_offset], ldvl)
00397                 ;
00398         i__1 = irows - 1;
00399         i__2 = irows - 1;
00400         template_lapack_lacpy("L", &i__1, &i__2, &b_ref(ilo + 1, ilo), ldb, &vl_ref(ilo + 1,
00401                  ilo), ldvl);
00402         i__1 = *lwork + 1 - iwrk;
00403         template_lapack_orgqr(&irows, &irows, &irows, &vl_ref(ilo, ilo), ldvl, &work[itau], 
00404                 &work[iwrk], &i__1, &ierr);
00405     }
00406 
00407 /*     Initialize VR */
00408 
00409     if (ilvr) {
00410         template_lapack_laset("Full", n, n, &c_b26, &c_b27, &vr[vr_offset], ldvr)
00411                 ;
00412     }
00413 
00414 /*     Reduce to generalized Hessenberg form   
00415        (Workspace: none needed) */
00416 
00417     if (ilv) {
00418 
00419 /*        Eigenvectors requested -- work on whole matrix. */
00420 
00421         template_lapack_gghrd(jobvl, jobvr, n, &ilo, &ihi, &a[a_offset], lda, &b[b_offset], 
00422                 ldb, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, &ierr);
00423     } else {
00424         template_lapack_gghrd("N", "N", &irows, &c__1, &irows, &a_ref(ilo, ilo), lda, &
00425                 b_ref(ilo, ilo), ldb, &vl[vl_offset], ldvl, &vr[vr_offset], 
00426                 ldvr, &ierr);
00427     }
00428 
00429 /*     Perform QZ algorithm (Compute eigenvalues, and optionally, the   
00430        Schur forms and Schur vectors)   
00431        (Workspace: need N) */
00432 
00433     iwrk = itau;
00434     if (ilv) {
00435         *(unsigned char *)chtemp = 'S';
00436     } else {
00437         *(unsigned char *)chtemp = 'E';
00438     }
00439     i__1 = *lwork + 1 - iwrk;
00440     template_lapack_hgeqz(chtemp, jobvl, jobvr, n, &ilo, &ihi, &a[a_offset], lda, &b[
00441             b_offset], ldb, &alphar[1], &alphai[1], &beta[1], &vl[vl_offset], 
00442             ldvl, &vr[vr_offset], ldvr, &work[iwrk], &i__1, &ierr);
00443     if (ierr != 0) {
00444         if (ierr > 0 && ierr <= *n) {
00445             *info = ierr;
00446         } else if (ierr > *n && ierr <= *n << 1) {
00447             *info = ierr - *n;
00448         } else {
00449             *info = *n + 1;
00450         }
00451         goto L110;
00452     }
00453 
00454 /*     Compute Eigenvectors   
00455        (Workspace: need 6*N) */
00456 
00457     if (ilv) {
00458         if (ilvl) {
00459             if (ilvr) {
00460                 *(unsigned char *)chtemp = 'B';
00461             } else {
00462                 *(unsigned char *)chtemp = 'L';
00463             }
00464         } else {
00465             *(unsigned char *)chtemp = 'R';
00466         }
00467         template_lapack_tgevc(chtemp, "B", ldumma, n, &a[a_offset], lda, &b[b_offset], ldb, 
00468                 &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, n, &in, &work[
00469                 iwrk], &ierr);
00470         if (ierr != 0) {
00471             *info = *n + 2;
00472             goto L110;
00473         }
00474 
00475 /*        Undo balancing on VL and VR and normalization   
00476           (Workspace: none needed) */
00477 
00478         if (ilvl) {
00479             template_lapack_ggbak("P", "L", n, &ilo, &ihi, &work[ileft], &work[iright], n, &
00480                     vl[vl_offset], ldvl, &ierr);
00481             i__1 = *n;
00482             for (jc = 1; jc <= i__1; ++jc) {
00483                 if (alphai[jc] < 0.) {
00484                     goto L50;
00485                 }
00486                 temp = 0.;
00487                 if (alphai[jc] == 0.) {
00488                     i__2 = *n;
00489                     for (jr = 1; jr <= i__2; ++jr) {
00490 /* Computing MAX */
00491                         d__2 = temp, d__3 = (d__1 = vl_ref(jr, jc), absMACRO(d__1))
00492                                 ;
00493                         temp = maxMACRO(d__2,d__3);
00494 /* L10: */
00495                     }
00496                 } else {
00497                     i__2 = *n;
00498                     for (jr = 1; jr <= i__2; ++jr) {
00499 /* Computing MAX */
00500                         d__3 = temp, d__4 = (d__1 = vl_ref(jr, jc), absMACRO(d__1))
00501                                  + (d__2 = vl_ref(jr, jc + 1), absMACRO(d__2));
00502                         temp = maxMACRO(d__3,d__4);
00503 /* L20: */
00504                     }
00505                 }
00506                 if (temp < smlnum) {
00507                     goto L50;
00508                 }
00509                 temp = 1. / temp;
00510                 if (alphai[jc] == 0.) {
00511                     i__2 = *n;
00512                     for (jr = 1; jr <= i__2; ++jr) {
00513                         vl_ref(jr, jc) = vl_ref(jr, jc) * temp;
00514 /* L30: */
00515                     }
00516                 } else {
00517                     i__2 = *n;
00518                     for (jr = 1; jr <= i__2; ++jr) {
00519                         vl_ref(jr, jc) = vl_ref(jr, jc) * temp;
00520                         vl_ref(jr, jc + 1) = vl_ref(jr, jc + 1) * temp;
00521 /* L40: */
00522                     }
00523                 }
00524 L50:
00525                 ;
00526             }
00527         }
00528         if (ilvr) {
00529             template_lapack_ggbak("P", "R", n, &ilo, &ihi, &work[ileft], &work[iright], n, &
00530                     vr[vr_offset], ldvr, &ierr);
00531             i__1 = *n;
00532             for (jc = 1; jc <= i__1; ++jc) {
00533                 if (alphai[jc] < 0.) {
00534                     goto L100;
00535                 }
00536                 temp = 0.;
00537                 if (alphai[jc] == 0.) {
00538                     i__2 = *n;
00539                     for (jr = 1; jr <= i__2; ++jr) {
00540 /* Computing MAX */
00541                         d__2 = temp, d__3 = (d__1 = vr_ref(jr, jc), absMACRO(d__1))
00542                                 ;
00543                         temp = maxMACRO(d__2,d__3);
00544 /* L60: */
00545                     }
00546                 } else {
00547                     i__2 = *n;
00548                     for (jr = 1; jr <= i__2; ++jr) {
00549 /* Computing MAX */
00550                         d__3 = temp, d__4 = (d__1 = vr_ref(jr, jc), absMACRO(d__1))
00551                                  + (d__2 = vr_ref(jr, jc + 1), absMACRO(d__2));
00552                         temp = maxMACRO(d__3,d__4);
00553 /* L70: */
00554                     }
00555                 }
00556                 if (temp < smlnum) {
00557                     goto L100;
00558                 }
00559                 temp = 1. / temp;
00560                 if (alphai[jc] == 0.) {
00561                     i__2 = *n;
00562                     for (jr = 1; jr <= i__2; ++jr) {
00563                         vr_ref(jr, jc) = vr_ref(jr, jc) * temp;
00564 /* L80: */
00565                     }
00566                 } else {
00567                     i__2 = *n;
00568                     for (jr = 1; jr <= i__2; ++jr) {
00569                         vr_ref(jr, jc) = vr_ref(jr, jc) * temp;
00570                         vr_ref(jr, jc + 1) = vr_ref(jr, jc + 1) * temp;
00571 /* L90: */
00572                     }
00573                 }
00574 L100:
00575                 ;
00576             }
00577         }
00578 
00579 /*        End of eigenvector calculation */
00580 
00581     }
00582 
00583 /*     Undo scaling if necessary */
00584 
00585     if (ilascl) {
00586         template_lapack_lascl("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphar[1], n, &
00587                 ierr);
00588         template_lapack_lascl("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphai[1], n, &
00589                 ierr);
00590     }
00591 
00592     if (ilbscl) {
00593         template_lapack_lascl("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, &beta[1], n, &
00594                 ierr);
00595     }
00596 
00597 L110:
00598 
00599     work[1] = (Treal) maxwrk;
00600 
00601     return 0;
00602 
00603 /*     End of DGGEV */
00604 
00605 } /* dggev_ */
00606 
00607 #undef vr_ref
00608 #undef vl_ref
00609 #undef b_ref
00610 #undef a_ref
00611 
00612 
00613 #endif

Generated on Mon Sep 17 14:30:40 2012 for ergo by  doxygen 1.4.7