/* ../SRC/zneupd.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    integer logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt, msapps, 
	    msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, 
	    mneupd, mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd;
} debug_;

#define debug_1 debug_

struct {
    integer nopx, nbx, nrorth, nitref, nrstrt;
    real tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd, 
	    tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, tcaupd, tcaup2, 
	    tcaitr, tceigh, tcgets, tcapps, tcconv, tmvopx, tmvbx, tgetv0, 
	    titref, trvec;
} timing_;

#define timing_1 timing_

/* Table of constant values */

static doublecomplex c_b1 = {1.,0.};
static doublecomplex c_b2 = {0.,0.};
static doublereal c_b5 = .66666666666666663;
static integer c__1 = 1;
static logical c_true = TRUE_;

/* \BeginDoc */

/* \Name: zneupd */

/* \Description: */
/*  This subroutine returns the converged approximations to eigenvalues */
/*  of A*z = lambda*B*z and (optionally): */

/*      (1) The corresponding approximate eigenvectors; */

/*      (2) An orthonormal basis for the associated approximate */
/*          invariant subspace; */

/*      (3) Both. */

/*  There is negligible additional cost to obtain eigenvectors.  An orthonormal */
/*  basis is always computed.  There is an additional storage cost of n*nev */
/*  if both are requested (in this case a separate array Z must be supplied). */

/*  The approximate eigenvalues and eigenvectors of  A*z = lambda*B*z */
/*  are derived from approximate eigenvalues and eigenvectors of */
/*  of the linear operator OP prescribed by the MODE selection in the */
/*  call to ZNAUPD.  ZNAUPD must be called before this routine is called. */
/*  These approximate eigenvalues and vectors are commonly called Ritz */
/*  values and Ritz vectors respectively.  They are referred to as such */
/*  in the comments that follow.   The computed orthonormal basis for the */
/*  invariant subspace corresponding to these Ritz values is referred to as a */
/*  Schur basis. */

/*  The definition of OP as well as other terms and the relation of computed */
/*  Ritz values and vectors of OP with respect to the given problem */
/*  A*z = lambda*B*z may be found in the header of ZNAUPD.  For a brief */
/*  description, see definitions of IPARAM(7), MODE and WHICH in the */
/*  documentation of ZNAUPD. */

/* \Usage: */
/*  call zneupd */
/*     ( RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, WORKEV, BMAT, */
/*       N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, */
/*       WORKL, LWORKL, RWORK, INFO ) */

/* \Arguments: */
/*  RVEC    LOGICAL  (INPUT) */
/*          Specifies whether a basis for the invariant subspace corresponding */
/*          to the converged Ritz value approximations for the eigenproblem */
/*          A*z = lambda*B*z is computed. */

/*             RVEC = .FALSE.     Compute Ritz values only. */

/*             RVEC = .TRUE.      Compute Ritz vectors or Schur vectors. */
/*                                See Remarks below. */

/*  HOWMNY  Character*1  (INPUT) */
/*          Specifies the form of the basis for the invariant subspace */
/*          corresponding to the converged Ritz values that is to be computed. */

/*          = 'A': Compute NEV Ritz vectors; */
/*          = 'P': Compute NEV Schur vectors; */
/*          = 'S': compute some of the Ritz vectors, specified */
/*                 by the logical array SELECT. */

/*  SELECT  Logical array of dimension NCV.  (INPUT) */
/*          If HOWMNY = 'S', SELECT specifies the Ritz vectors to be */
/*          computed. To select the  Ritz vector corresponding to a */
/*          Ritz value D(j), SELECT(j) must be set to .TRUE.. */
/*          If HOWMNY = 'A' or 'P', SELECT need not be initialized */
/*          but it is used as internal workspace. */

/*  D       Complex*16 array of dimension NEV+1.  (OUTPUT) */
/*          On exit, D contains the  Ritz  approximations */
/*          to the eigenvalues lambda for A*z = lambda*B*z. */

/*  Z       Complex*16 N by NEV array    (OUTPUT) */
/*          On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of */
/*          Z represents approximate eigenvectors (Ritz vectors) corresponding */
/*          to the NCONV=IPARAM(5) Ritz values for eigensystem */
/*          A*z = lambda*B*z. */

/*          If RVEC = .FALSE. or HOWMNY = 'P', then Z is NOT REFERENCED. */

/*          NOTE: If if RVEC = .TRUE. and a Schur basis is not required, */
/*          the array Z may be set equal to first NEV+1 columns of the Arnoldi */
/*          basis array V computed by ZNAUPD.  In this case the Arnoldi basis */
/*          will be destroyed and overwritten with the eigenvector basis. */

/*  LDZ     Integer.  (INPUT) */
/*          The leading dimension of the array Z.  If Ritz vectors are */
/*          desired, then  LDZ .ge.  max( 1, N ) is required. */
/*          In any case,  LDZ .ge. 1 is required. */

/*  SIGMA   Complex*16  (INPUT) */
/*          If IPARAM(7) = 3 then SIGMA represents the shift. */
/*          Not referenced if IPARAM(7) = 1 or 2. */

/*  WORKEV  Complex*16 work array of dimension 2*NCV.  (WORKSPACE) */

/*  **** The remaining arguments MUST be the same as for the   **** */
/*  **** call to ZNAUPD that was just completed.               **** */

/*  NOTE: The remaining arguments */

/*           BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, */
/*           WORKD, WORKL, LWORKL, RWORK, INFO */

/*         must be passed directly to ZNEUPD following the last call */
/*         to ZNAUPD.  These arguments MUST NOT BE MODIFIED between */
/*         the the last call to ZNAUPD and the call to ZNEUPD. */

/*  Three of these parameters (V, WORKL and INFO) are also output parameters: */

/*  V       Complex*16 N by NCV array.  (INPUT/OUTPUT) */

/*          Upon INPUT: the NCV columns of V contain the Arnoldi basis */
/*                      vectors for OP as constructed by ZNAUPD . */

/*          Upon OUTPUT: If RVEC = .TRUE. the first NCONV=IPARAM(5) columns */
/*                       contain approximate Schur vectors that span the */
/*                       desired invariant subspace. */

/*          NOTE: If the array Z has been set equal to first NEV+1 columns */
/*          of the array V and RVEC=.TRUE. and HOWMNY= 'A', then the */
/*          Arnoldi basis held by V has been overwritten by the desired */
/*          Ritz vectors.  If a separate array Z has been passed then */
/*          the first NCONV=IPARAM(5) columns of V will contain approximate */
/*          Schur vectors that span the desired invariant subspace. */

/*  WORKL   Double precision work array of length LWORKL.  (OUTPUT/WORKSPACE) */
/*          WORKL(1:ncv*ncv+2*ncv) contains information obtained in */
/*          znaupd.  They are not changed by zneupd. */
/*          WORKL(ncv*ncv+2*ncv+1:3*ncv*ncv+4*ncv) holds the */
/*          untransformed Ritz values, the untransformed error estimates of */
/*          the Ritz values, the upper triangular matrix for H, and the */
/*          associated matrix representation of the invariant subspace for H. */

/*          Note: IPNTR(9:13) contains the pointer into WORKL for addresses */
/*          of the above information computed by zneupd. */
/*          ------------------------------------------------------------- */
/*          IPNTR(9):  pointer to the NCV RITZ values of the */
/*                     original system. */
/*          IPNTR(10): Not used */
/*          IPNTR(11): pointer to the NCV corresponding error estimates. */
/*          IPNTR(12): pointer to the NCV by NCV upper triangular */
/*                     Schur matrix for H. */
/*          IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors */
/*                     of the upper Hessenberg matrix H. Only referenced by */
/*                     zneupd if RVEC = .TRUE. See Remark 2 below. */
/*          ------------------------------------------------------------- */

/*  INFO    Integer.  (OUTPUT) */
/*          Error flag on output. */
/*          =  0: Normal exit. */

/*          =  1: The Schur form computed by LAPACK routine csheqr */
/*                could not be reordered by LAPACK routine ztrsen. */
/*                Re-enter subroutine zneupd with IPARAM(5)=NCV and */
/*                increase the size of the array D to have */
/*                dimension at least dimension NCV and allocate at least NCV */
/*                columns for Z. NOTE: Not necessary if Z and V share */
/*                the same space. Please notify the authors if this error */
/*                occurs. */

/*          = -1: N must be positive. */
/*          = -2: NEV must be positive. */
/*          = -3: NCV-NEV >= 2 and less than or equal to N. */
/*          = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI' */
/*          = -6: BMAT must be one of 'I' or 'G'. */
/*          = -7: Length of private work WORKL array is not sufficient. */
/*          = -8: Error return from LAPACK eigenvalue calculation. */
/*                This should never happened. */
/*          = -9: Error return from calculation of eigenvectors. */
/*                Informational error from LAPACK routine ztrevc. */
/*          = -10: IPARAM(7) must be 1,2,3 */
/*          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible. */
/*          = -12: HOWMNY = 'S' not yet implemented */
/*          = -13: HOWMNY must be one of 'A' or 'P' if RVEC = .true. */
/*          = -14: ZNAUPD did not find any eigenvalues to sufficient */
/*                 accuracy. */

/* \BeginLib */

/* \References: */
/*  1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
/*     a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
/*     pp 357-385. */
/*  2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
/*     Restarted Arnoldi Iteration", Rice University Technical Report */
/*     TR95-13, Department of Computational and Applied Mathematics. */
/*  3. B. Nour-Omid, B. N. Parlett, T. Ericsson and P. S. Jensen, */
/*     "How to Implement the Spectral Transformation", Math Comp., */
/*     Vol. 48, No. 178, April, 1987 pp. 664-673. */

/* \Routines called: */
/*     ivout   ARPACK utility routine that prints integers. */
/*     zmout   ARPACK utility routine that prints matrices */
/*     zvout   ARPACK utility routine that prints vectors. */
/*     zgeqr2  LAPACK routine that computes the QR factorization of */
/*             a matrix. */
/*     zlacpy  LAPACK matrix copy routine. */
/*     zlahqr  LAPACK routine that computes the Schur form of a */
/*             upper Hessenberg matrix. */
/*     zlaset  LAPACK matrix initialization routine. */
/*     ztrevc  LAPACK routine to compute the eigenvectors of a matrix */
/*             in upper triangular form. */
/*     ztrsen  LAPACK routine that re-orders the Schur form. */
/*     zunm2r  LAPACK routine that applies an orthogonal matrix in */
/*             factored form. */
/*     dlamch  LAPACK routine that determines machine constants. */
/*     ztrmm   Level 3 BLAS matrix times an upper triangular matrix. */
/*     zgeru   Level 2 BLAS rank one update to a matrix. */
/*     zcopy   Level 1 BLAS that copies one vector to another . */
/*     zscal   Level 1 BLAS that scales a vector. */
/*     zdscal  Level 1 BLAS that scales a complex vector by a real number. */
/*     dznrm2  Level 1 BLAS that computes the norm of a complex vector. */

/* \Remarks */

/*  1. Currently only HOWMNY = 'A' and 'P' are implemented. */

/*  2. Schur vectors are an orthogonal representation for the basis of */
/*     Ritz vectors. Thus, their numerical properties are often superior. */
/*     If RVEC = .true. then the relationship */
/*             A * V(:,1:IPARAM(5)) = V(:,1:IPARAM(5)) * T, and */
/*     V(:,1:IPARAM(5))' * V(:,1:IPARAM(5)) = I are approximately satisfied. */
/*     Here T is the leading submatrix of order IPARAM(5) of the */
/*     upper triangular matrix stored workl(ipntr(12)). */

/* \Authors */
/*     Danny Sorensen               Phuong Vu */
/*     Richard Lehoucq              CRPC / Rice University */
/*     Chao Yang                    Houston, Texas */
/*     Dept. of Computational & */
/*     Applied Mathematics */
/*     Rice University */
/*     Houston, Texas */

/* \SCCS Information: @(#) */
/* FILE: neupd.F   SID: 2.4   DATE OF SID: 7/31/96   RELEASE: 2 */

/* \EndLib */

/* ----------------------------------------------------------------------- */
/* Subroutine */ int zneupd_(logical *rvec, char *howmny, logical *select, 
	doublecomplex *d__, doublecomplex *z__, integer *ldz, doublecomplex *
	sigma, doublecomplex *workev, char *bmat, integer *n, char *which, 
	integer *nev, doublereal *tol, doublecomplex *resid, integer *ncv, 
	doublecomplex *v, integer *ldv, integer *iparam, integer *ipntr, 
	doublecomplex *workd, doublecomplex *workl, integer *lworkl, 
	doublereal *rwork, integer *info, ftnlen howmny_len, ftnlen bmat_len, 
	ftnlen which_len)
{
    /* System generated locals */
    integer v_dim1, v_offset, z_dim1, z_offset, i__1, i__2;
    doublereal d__1, d__2, d__3, d__4;
    doublecomplex z__1, z__2;

    /* Builtin functions */
    double pow_dd(doublereal *, doublereal *);
    integer s_cmp(char *, char *, ftnlen, ftnlen);
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
    double d_imag(doublecomplex *);
    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);

    /* Local variables */
    static integer j, k, ih, iq;
    static doublecomplex vl[1];
    static integer wr, ibd, ldh, ldq;
    static doublereal sep;
    static integer irz, mode;
    static doublereal eps23;
    static integer ierr;
    static doublecomplex temp;
    static integer iwev;
    static char type__[6];
    static integer ritz, iheig, ihbds;
    static doublereal conds;
    static logical reord;
    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 
	    doublecomplex *, integer *);
    static integer nconv;
    static doublereal thres, rtemp;
    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 
	    doublecomplex *, integer *, doublecomplex *, integer *);
    static doublecomplex rnorm;
    extern /* Subroutine */ int zgeru_(integer *, integer *, doublecomplex *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, 
	    doublecomplex *, integer *), dvout_(integer *, integer *, 
	    doublereal *, integer *, char *, ftnlen), zcopy_(integer *, 
	    doublecomplex *, integer *, doublecomplex *, integer *), ivout_(
	    integer *, integer *, integer *, integer *, char *, ftnlen), 
	    ztrmm_(char *, char *, char *, char *, integer *, integer *, 
	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
	    integer *, ftnlen, ftnlen, ftnlen, ftnlen), zmout_(integer *, 
	    integer *, integer *, doublecomplex *, integer *, integer *, char 
	    *, ftnlen), zvout_(integer *, integer *, doublecomplex *, integer 
	    *, char *, ftnlen);
    extern doublereal dlapy2_(doublereal *, doublereal *);
    extern /* Subroutine */ int zgeqr2_(integer *, integer *, doublecomplex *,
	     integer *, doublecomplex *, doublecomplex *, integer *);
    extern doublereal dznrm2_(integer *, doublecomplex *, integer *), dlamch_(
	    char *, ftnlen);
    extern /* Subroutine */ int zunm2r_(char *, char *, integer *, integer *, 
	    integer *, doublecomplex *, integer *, doublecomplex *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, ftnlen, 
	    ftnlen);
    static integer bounds, invsub, iuptri, msglvl, ktrord, outncv;
    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, ftnlen), 
	    zlahqr_(logical *, logical *, integer *, integer *, integer *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, integer *,
	     doublecomplex *, integer *, integer *), zlaset_(char *, integer *
	    , integer *, doublecomplex *, doublecomplex *, doublecomplex *, 
	    integer *, ftnlen), ztrsen_(char *, char *, logical *, integer *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, 
	    doublecomplex *, integer *, doublereal *, doublereal *, 
	    doublecomplex *, integer *, integer *, ftnlen, ftnlen), ztrevc_(
	    char *, char *, logical *, integer *, doublecomplex *, integer *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, integer *,
	     integer *, doublecomplex *, doublereal *, integer *, ftnlen, 
	    ftnlen), zdscal_(integer *, doublereal *, doublecomplex *, 
	    integer *);


/*     %----------------------------------------------------% */
/*     | Include files for debugging and timing information | */
/*     %----------------------------------------------------% */


/* \SCCS Information: @(#) */
/* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 */

/*     %---------------------------------% */
/*     | See debug.doc for documentation | */
/*     %---------------------------------% */

/*     %------------------% */
/*     | Scalar Arguments | */
/*     %------------------% */

/*     %--------------------------------% */
/*     | See stat.doc for documentation | */
/*     %--------------------------------% */

/* \SCCS Information: @(#) */
/* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 */



/*     %-----------------% */
/*     | Array Arguments | */
/*     %-----------------% */


/*     %------------% */
/*     | Parameters | */
/*     %------------% */


/*     %---------------% */
/*     | Local Scalars | */
/*     %---------------% */


/*     %----------------------% */
/*     | External Subroutines | */
/*     %----------------------% */


/*     %--------------------% */
/*     | External Functions | */
/*     %--------------------% */



/*     %-----------------------% */
/*     | Executable Statements | */
/*     %-----------------------% */

/*     %------------------------% */
/*     | Set default parameters | */
/*     %------------------------% */

    /* Parameter adjustments */
    --workd;
    --resid;
    z_dim1 = *ldz;
    z_offset = 1 + z_dim1;
    z__ -= z_offset;
    --d__;
    --rwork;
    --workev;
    --select;
    v_dim1 = *ldv;
    v_offset = 1 + v_dim1;
    v -= v_offset;
    --iparam;
    --ipntr;
    --workl;

    /* Function Body */
    msglvl = debug_1.mceupd;
    mode = iparam[7];
    nconv = iparam[5];
    *info = 0;


/*     %---------------------------------% */
/*     | Get machine dependent constant. | */
/*     %---------------------------------% */

    eps23 = dlamch_("Epsilon-Machine", (ftnlen)15);
    eps23 = pow_dd(&eps23, &c_b5);

/*     %-------------------------------% */
/*     | Quick return                  | */
/*     | Check for incompatible input  | */
/*     %-------------------------------% */

    ierr = 0;

    if (nconv <= 0) {
	ierr = -14;
    } else if (*n <= 0) {
	ierr = -1;
    } else if (*nev <= 0) {
	ierr = -2;
    } else if (*ncv <= *nev + 1 || *ncv > *n) {
	ierr = -3;
    } else if (s_cmp(which, "LM", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which, 
	    "SM", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which, "LR", (ftnlen)2, 
	    (ftnlen)2) != 0 && s_cmp(which, "SR", (ftnlen)2, (ftnlen)2) != 0 
	    && s_cmp(which, "LI", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which, 
	    "SI", (ftnlen)2, (ftnlen)2) != 0) {
	ierr = -5;
    } else if (*(unsigned char *)bmat != 'I' && *(unsigned char *)bmat != 'G')
	     {
	ierr = -6;
    } else /* if(complicated condition) */ {
/* Computing 2nd power */
	i__1 = *ncv;
	if (*lworkl < i__1 * i__1 * 3 + (*ncv << 2)) {
	    ierr = -7;
	} else if (*(unsigned char *)howmny != 'A' && *(unsigned char *)
		howmny != 'P' && *(unsigned char *)howmny != 'S' && *rvec) {
	    ierr = -13;
	} else if (*(unsigned char *)howmny == 'S') {
	    ierr = -12;
	}
    }

    if (mode == 1 || mode == 2) {
	s_copy(type__, "REGULR", (ftnlen)6, (ftnlen)6);
    } else if (mode == 3) {
	s_copy(type__, "SHIFTI", (ftnlen)6, (ftnlen)6);
    } else {
	ierr = -10;
    }
    if (mode == 1 && *(unsigned char *)bmat == 'G') {
	ierr = -11;
    }

/*     %------------% */
/*     | Error Exit | */
/*     %------------% */

    if (ierr != 0) {
	*info = ierr;
	goto L9000;
    }

/*     %--------------------------------------------------------% */
/*     | Pointer into WORKL for address of H, RITZ, WORKEV, Q   | */
/*     | etc... and the remaining workspace.                    | */
/*     | Also update pointer to be used on output.              | */
/*     | Memory is laid out as follows:                         | */
/*     | workl(1:ncv*ncv) := generated Hessenberg matrix        | */
/*     | workl(ncv*ncv+1:ncv*ncv+ncv) := ritz values            | */
/*     | workl(ncv*ncv+ncv+1:ncv*ncv+2*ncv) := error bounds     | */
/*     %--------------------------------------------------------% */

/*     %-----------------------------------------------------------% */
/*     | The following is used and set by ZNEUPD.                  | */
/*     | workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := The untransformed | */
/*     |                                      Ritz values.         | */
/*     | workl(ncv*ncv+3*ncv+1:ncv*ncv+4*ncv) := The untransformed | */
/*     |                                      error bounds of      | */
/*     |                                      the Ritz values      | */
/*     | workl(ncv*ncv+4*ncv+1:2*ncv*ncv+4*ncv) := Holds the upper | */
/*     |                                      triangular matrix    | */
/*     |                                      for H.               | */
/*     | workl(2*ncv*ncv+4*ncv+1: 3*ncv*ncv+4*ncv) := Holds the    | */
/*     |                                      associated matrix    | */
/*     |                                      representation of    | */
/*     |                                      the invariant        | */
/*     |                                      subspace for H.      | */
/*     | GRAND total of NCV * ( 3 * NCV + 4 ) locations.           | */
/*     %-----------------------------------------------------------% */

    ih = ipntr[5];
    ritz = ipntr[6];
    iq = ipntr[7];
    bounds = ipntr[8];
    ldh = *ncv;
    ldq = *ncv;
    iheig = bounds + ldh;
    ihbds = iheig + ldh;
    iuptri = ihbds + ldh;
    invsub = iuptri + ldh * *ncv;
    ipntr[9] = iheig;
    ipntr[11] = ihbds;
    ipntr[12] = iuptri;
    ipntr[13] = invsub;
    wr = 1;
    iwev = wr + *ncv;

/*     %-----------------------------------------% */
/*     | irz points to the Ritz values computed  | */
/*     |     by _neigh before exiting _naup2.    | */
/*     | ibd points to the Ritz estimates        | */
/*     |     computed by _neigh before exiting   | */
/*     |     _naup2.                             | */
/*     %-----------------------------------------% */

    irz = ipntr[14] + *ncv * *ncv;
    ibd = irz + *ncv;

/*     %------------------------------------% */
/*     | RNORM is B-norm of the RESID(1:N). | */
/*     %------------------------------------% */

    i__1 = ih + 2;
    rnorm.r = workl[i__1].r, rnorm.i = workl[i__1].i;
    i__1 = ih + 2;
    workl[i__1].r = 0., workl[i__1].i = 0.;

    if (*rvec) {

/*        %-------------------------------------------% */
/*        | Get converged Ritz value on the boundary. | */
/*        | Note: converged Ritz values have been     | */
/*        | placed in the first NCONV locations in    | */
/*        | workl(ritz).  They have been sorted       | */
/*        | (in _naup2) according to the WHICH        | */
/*        | selection criterion                       | */
/*        %-------------------------------------------% */

	if (s_cmp(which, "LM", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(which, 
		"SM", (ftnlen)2, (ftnlen)2) == 0) {
	    i__1 = ritz;
	    d__1 = workl[i__1].r;
	    d__2 = d_imag(&workl[ritz]);
	    thres = dlapy2_(&d__1, &d__2);
	} else if (s_cmp(which, "LR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
		which, "SR", (ftnlen)2, (ftnlen)2) == 0) {
	    i__1 = ritz;
	    thres = workl[i__1].r;
	} else if (s_cmp(which, "LI", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
		which, "SI", (ftnlen)2, (ftnlen)2) == 0) {
	    thres = d_imag(&workl[ritz]);
	}
	if (msglvl > 2) {
	    dvout_(&debug_1.logfil, &c__1, &thres, &debug_1.ndigit, "_neupd:"
		    " Threshold eigenvalue used for re-ordering", (ftnlen)49);
	}

/*        %---------------------------------------------------------% */
/*        | Check to see if all converged Ritz values appear at the | */
/*        | at the top of the upper triangular matrix computed by   | */
/*        | _neigh in _naup2.  This is done in the following way:   | */
/*        |                                                         | */
/*        | 1) For each Ritz value from _neigh, compare it with the | */
/*        |    threshold Ritz value computed above to determine     | */
/*        |    whether it is a wanted one.                          | */
/*        |                                                         | */
/*        | 2) If it is wanted, then check the corresponding Ritz   | */
/*        |    estimate to see if it has converged.  If it has, set | */
/*        |    correponding entry in the logical array SELECT to    | */
/*        |    .TRUE..                                              | */
/*        |                                                         | */
/*        | If SELECT(j) = .TRUE. and j > NCONV, then there is a    | */
/*        | converged Ritz value that does not appear at the top of | */
/*        | the upper triangular matrix computed by _neigh in       | */
/*        | _naup2.  Reordering is needed.                          | */
/*        %---------------------------------------------------------% */

	reord = FALSE_;
	ktrord = 0;
	i__1 = *ncv - 1;
	for (j = 0; j <= i__1; ++j) {
	    select[j + 1] = FALSE_;
	    if (s_cmp(which, "LM", (ftnlen)2, (ftnlen)2) == 0) {
		i__2 = irz + j;
		d__1 = workl[i__2].r;
		d__2 = d_imag(&workl[irz + j]);
		if (dlapy2_(&d__1, &d__2) >= thres) {
/* Computing MAX */
		    i__2 = irz + j - 1;
		    d__3 = workl[i__2].r;
		    d__4 = d_imag(&workl[irz + j - 1]);
		    d__1 = eps23, d__2 = dlapy2_(&d__3, &d__4);
		    rtemp = max(d__1,d__2);
		    i__2 = ibd + j;
		    d__1 = workl[i__2].r;
		    d__2 = d_imag(&workl[ibd + j]);
		    if (dlapy2_(&d__1, &d__2) <= *tol * rtemp) {
			select[j + 1] = TRUE_;
		    }
		}
	    } else if (s_cmp(which, "SM", (ftnlen)2, (ftnlen)2) == 0) {
		i__2 = irz + j;
		d__1 = workl[i__2].r;
		d__2 = d_imag(&workl[irz + j]);
		if (dlapy2_(&d__1, &d__2) <= thres) {
/* Computing MAX */
		    i__2 = irz + j - 1;
		    d__3 = workl[i__2].r;
		    d__4 = d_imag(&workl[irz + j - 1]);
		    d__1 = eps23, d__2 = dlapy2_(&d__3, &d__4);
		    rtemp = max(d__1,d__2);
		    i__2 = ibd + j;
		    d__1 = workl[i__2].r;
		    d__2 = d_imag(&workl[ibd + j]);
		    if (dlapy2_(&d__1, &d__2) <= *tol * rtemp) {
			select[j + 1] = TRUE_;
		    }
		}
	    } else if (s_cmp(which, "LR", (ftnlen)2, (ftnlen)2) == 0) {
		i__2 = irz + j;
		if (workl[i__2].r >= thres) {
/* Computing MAX */
		    i__2 = irz + j - 1;
		    d__3 = workl[i__2].r;
		    d__4 = d_imag(&workl[irz + j - 1]);
		    d__1 = eps23, d__2 = dlapy2_(&d__3, &d__4);
		    rtemp = max(d__1,d__2);
		    i__2 = ibd + j;
		    d__1 = workl[i__2].r;
		    d__2 = d_imag(&workl[ibd + j]);
		    if (dlapy2_(&d__1, &d__2) <= *tol * rtemp) {
			select[j + 1] = TRUE_;
		    }
		}
	    } else if (s_cmp(which, "SR", (ftnlen)2, (ftnlen)2) == 0) {
		i__2 = irz + j;
		if (workl[i__2].r <= thres) {
/* Computing MAX */
		    i__2 = irz + j - 1;
		    d__3 = workl[i__2].r;
		    d__4 = d_imag(&workl[irz + j - 1]);
		    d__1 = eps23, d__2 = dlapy2_(&d__3, &d__4);
		    rtemp = max(d__1,d__2);
		    i__2 = ibd + j;
		    d__1 = workl[i__2].r;
		    d__2 = d_imag(&workl[ibd + j]);
		    if (dlapy2_(&d__1, &d__2) <= *tol * rtemp) {
			select[j + 1] = TRUE_;
		    }
		}
	    } else if (s_cmp(which, "LI", (ftnlen)2, (ftnlen)2) == 0) {
		if (d_imag(&workl[irz + j]) >= thres) {
/* Computing MAX */
		    i__2 = irz + j - 1;
		    d__3 = workl[i__2].r;
		    d__4 = d_imag(&workl[irz + j - 1]);
		    d__1 = eps23, d__2 = dlapy2_(&d__3, &d__4);
		    rtemp = max(d__1,d__2);
		    i__2 = ibd + j;
		    d__1 = workl[i__2].r;
		    d__2 = d_imag(&workl[ibd + j]);
		    if (dlapy2_(&d__1, &d__2) <= *tol * rtemp) {
			select[j + 1] = TRUE_;
		    }
		}
	    } else if (s_cmp(which, "SI", (ftnlen)2, (ftnlen)2) == 0) {
		if (d_imag(&workl[irz + j]) <= thres) {
/* Computing MAX */
		    i__2 = irz + j - 1;
		    d__3 = workl[i__2].r;
		    d__4 = d_imag(&workl[irz + j - 1]);
		    d__1 = eps23, d__2 = dlapy2_(&d__3, &d__4);
		    rtemp = max(d__1,d__2);
		    i__2 = ibd + j;
		    d__1 = workl[i__2].r;
		    d__2 = d_imag(&workl[ibd + j]);
		    if (dlapy2_(&d__1, &d__2) <= *tol * rtemp) {
			select[j + 1] = TRUE_;
		    }
		}
	    }
	    if (j + 1 > nconv) {
		reord = select[j + 1] || reord;
	    }
	    if (select[j + 1]) {
		++ktrord;
	    }
/* L10: */
	}

	if (msglvl > 2) {
	    ivout_(&debug_1.logfil, &c__1, &ktrord, &debug_1.ndigit, "_neupd"
		    ": Number of specified eigenvalues", (ftnlen)39);
	    ivout_(&debug_1.logfil, &c__1, &nconv, &debug_1.ndigit, "_neupd:"
		    " Number of \"converged\" eigenvalues", (ftnlen)41);
	}

/*        if (ktrord .gt. nconv) then */

/*           %-----------------------------------% */
/*           | More than NCONV Ritz values have  | */
/*           | "converged", and they all satisfy | */
/*           | the WHICH selection criterion.    | */
/*           %-----------------------------------% */

/*           iparam(6) = ktrord */

/*        end if */

/*        %-------------------------------------------------------% */
/*        | Call LAPACK routine zlahqr to compute the Schur form  | */
/*        | of the upper Hessenberg matrix returned by ZNAUPD.    | */
/*        | Make a copy of the upper Hessenberg matrix.           | */
/*        | Initialize the Schur vector matrix Q to the identity. | */
/*        %-------------------------------------------------------% */

	i__1 = ldh * *ncv;
	zcopy_(&i__1, &workl[ih], &c__1, &workl[iuptri], &c__1);
	zlaset_("All", ncv, ncv, &c_b2, &c_b1, &workl[invsub], &ldq, (ftnlen)
		3);
	zlahqr_(&c_true, &c_true, ncv, &c__1, ncv, &workl[iuptri], &ldh, &
		workl[iheig], &c__1, ncv, &workl[invsub], &ldq, &ierr);
	zcopy_(ncv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds], &c__1);

	if (ierr != 0) {
	    *info = -8;
	    goto L9000;
	}

	if (msglvl > 1) {
	    zvout_(&debug_1.logfil, ncv, &workl[iheig], &debug_1.ndigit, 
		    "_neupd: Eigenvalues of H", (ftnlen)24);
	    zvout_(&debug_1.logfil, ncv, &workl[ihbds], &debug_1.ndigit, 
		    "_neupd: Last row of the Schur vector matrix", (ftnlen)43)
		    ;
	    if (msglvl > 3) {
		zmout_(&debug_1.logfil, ncv, ncv, &workl[iuptri], &ldh, &
			debug_1.ndigit, "_neupd: The upper triangular matrix "
			, (ftnlen)36);
	    }
	}

	if (reord) {

/*           %-----------------------------------------------% */
/*           | Reorder the computed upper triangular matrix. | */
/*           %-----------------------------------------------% */

	    ztrsen_("None", "V", &select[1], ncv, &workl[iuptri], &ldh, &
		    workl[invsub], &ldq, &workl[iheig], &nconv, &conds, &sep, 
		    &workev[1], ncv, &ierr, (ftnlen)4, (ftnlen)1);

	    if (ierr == 1) {
		*info = 1;
		goto L9000;
	    }

	    if (msglvl > 2) {
		zvout_(&debug_1.logfil, ncv, &workl[iheig], &debug_1.ndigit, 
			"_neupd: Eigenvalues of H--reordered", (ftnlen)35);
		if (msglvl > 3) {
		    zmout_(&debug_1.logfil, ncv, ncv, &workl[iuptri], &ldq, &
			    debug_1.ndigit, "_neupd: Triangular matrix after"
			    " re-ordering", (ftnlen)43);
		}
	    }

	}

/*        %---------------------------------------------% */
/*        | Copy the last row of the Schur basis matrix | */
/*        | to workl(ihbds).  This vector will be used  | */
/*        | to compute the Ritz estimates of converged  | */
/*        | Ritz values.                                | */
/*        %---------------------------------------------% */

	zcopy_(ncv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds], &c__1);

/*        %--------------------------------------------% */
/*        | Place the computed eigenvalues of H into D | */
/*        | if a spectral transformation was not used. | */
/*        %--------------------------------------------% */

	if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) == 0) {
	    zcopy_(&nconv, &workl[iheig], &c__1, &d__[1], &c__1);
	}

/*        %----------------------------------------------------------% */
/*        | Compute the QR factorization of the matrix representing  | */
/*        | the wanted invariant subspace located in the first NCONV | */
/*        | columns of workl(invsub,ldq).                            | */
/*        %----------------------------------------------------------% */

	zgeqr2_(ncv, &nconv, &workl[invsub], &ldq, &workev[1], &workev[*ncv + 
		1], &ierr);

/*        %--------------------------------------------------------% */
/*        | * Postmultiply V by Q using zunm2r.                    | */
/*        | * Copy the first NCONV columns of VQ into Z.           | */
/*        | * Postmultiply Z by R.                                 | */
/*        | The N by NCONV matrix Z is now a matrix representation | */
/*        | of the approximate invariant subspace associated with  | */
/*        | the Ritz values in workl(iheig). The first NCONV       | */
/*        | columns of V are now approximate Schur vectors         | */
/*        | associated with the upper triangular matrix of order   | */
/*        | NCONV in workl(iuptri).                                | */
/*        %--------------------------------------------------------% */

	zunm2r_("Right", "Notranspose", n, ncv, &nconv, &workl[invsub], &ldq, 
		&workev[1], &v[v_offset], ldv, &workd[*n + 1], &ierr, (ftnlen)
		5, (ftnlen)11);
	zlacpy_("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz, (
		ftnlen)3);

	i__1 = nconv;
	for (j = 1; j <= i__1; ++j) {

/*           %---------------------------------------------------% */
/*           | Perform both a column and row scaling if the      | */
/*           | diagonal element of workl(invsub,ldq) is negative | */
/*           | I'm lazy and don't take advantage of the upper    | */
/*           | triangular form of workl(iuptri,ldq).             | */
/*           | Note that since Q is orthogonal, R is a diagonal  | */
/*           | matrix consisting of plus or minus ones.          | */
/*           %---------------------------------------------------% */

	    i__2 = invsub + (j - 1) * ldq + j - 1;
	    if (workl[i__2].r < 0.) {
		z__1.r = -1., z__1.i = -0.;
		zscal_(&nconv, &z__1, &workl[iuptri + j - 1], &ldq);
		z__1.r = -1., z__1.i = -0.;
		zscal_(&nconv, &z__1, &workl[iuptri + (j - 1) * ldq], &c__1);
	    }

/* L20: */
	}

	if (*(unsigned char *)howmny == 'A') {

/*           %--------------------------------------------% */
/*           | Compute the NCONV wanted eigenvectors of T | */
/*           | located in workl(iuptri,ldq).              | */
/*           %--------------------------------------------% */

	    i__1 = *ncv;
	    for (j = 1; j <= i__1; ++j) {
		if (j <= nconv) {
		    select[j] = TRUE_;
		} else {
		    select[j] = FALSE_;
		}
/* L30: */
	    }

	    ztrevc_("Right", "Select", &select[1], ncv, &workl[iuptri], &ldq, 
		    vl, &c__1, &workl[invsub], &ldq, ncv, &outncv, &workev[1],
		     &rwork[1], &ierr, (ftnlen)5, (ftnlen)6);

	    if (ierr != 0) {
		*info = -9;
		goto L9000;
	    }

/*           %------------------------------------------------% */
/*           | Scale the returning eigenvectors so that their | */
/*           | Euclidean norms are all one. LAPACK subroutine | */
/*           | ztrevc returns each eigenvector normalized so  | */
/*           | that the element of largest magnitude has      | */
/*           | magnitude 1.                                   | */
/*           %------------------------------------------------% */

	    i__1 = nconv;
	    for (j = 1; j <= i__1; ++j) {
		rtemp = dznrm2_(ncv, &workl[invsub + (j - 1) * ldq], &c__1);
		rtemp = 1. / rtemp;
		zdscal_(ncv, &rtemp, &workl[invsub + (j - 1) * ldq], &c__1);

/*                 %------------------------------------------% */
/*                 | Ritz estimates can be obtained by taking | */
/*                 | the inner product of the last row of the | */
/*                 | Schur basis of H with eigenvectors of T. | */
/*                 | Note that the eigenvector matrix of T is | */
/*                 | upper triangular, thus the length of the | */
/*                 | inner product can be set to j.           | */
/*                 %------------------------------------------% */

		i__2 = j;
		zdotc_(&z__1, &j, &workl[ihbds], &c__1, &workl[invsub + (j - 
			1) * ldq], &c__1);
		workev[i__2].r = z__1.r, workev[i__2].i = z__1.i;
/* L40: */
	    }

	    if (msglvl > 2) {
		zcopy_(&nconv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds],
			 &c__1);
		zvout_(&debug_1.logfil, &nconv, &workl[ihbds], &
			debug_1.ndigit, "_neupd: Last row of the eigenvector"
			" matrix for T", (ftnlen)48);
		if (msglvl > 3) {
		    zmout_(&debug_1.logfil, ncv, ncv, &workl[invsub], &ldq, &
			    debug_1.ndigit, "_neupd: The eigenvector matrix "
			    "for T", (ftnlen)36);
		}
	    }

/*           %---------------------------------------% */
/*           | Copy Ritz estimates into workl(ihbds) | */
/*           %---------------------------------------% */

	    zcopy_(&nconv, &workev[1], &c__1, &workl[ihbds], &c__1);

/*           %----------------------------------------------% */
/*           | The eigenvector matrix Q of T is triangular. | */
/*           | Form Z*Q.                                    | */
/*           %----------------------------------------------% */

	    ztrmm_("Right", "Upper", "No transpose", "Non-unit", n, &nconv, &
		    c_b1, &workl[invsub], &ldq, &z__[z_offset], ldz, (ftnlen)
		    5, (ftnlen)5, (ftnlen)12, (ftnlen)8);

	}

    } else {

/*        %--------------------------------------------------% */
/*        | An approximate invariant subspace is not needed. | */
/*        | Place the Ritz values computed ZNAUPD into D.    | */
/*        %--------------------------------------------------% */

	zcopy_(&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
	zcopy_(&nconv, &workl[ritz], &c__1, &workl[iheig], &c__1);
	zcopy_(&nconv, &workl[bounds], &c__1, &workl[ihbds], &c__1);

    }

/*     %------------------------------------------------% */
/*     | Transform the Ritz values and possibly vectors | */
/*     | and corresponding error bounds of OP to those  | */
/*     | of A*x = lambda*B*x.                           | */
/*     %------------------------------------------------% */

    if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) == 0) {

	if (*rvec) {
	    zscal_(ncv, &rnorm, &workl[ihbds], &c__1);
	}

    } else {

/*        %---------------------------------------% */
/*        |   A spectral transformation was used. | */
/*        | * Determine the Ritz estimates of the | */
/*        |   Ritz values in the original system. | */
/*        %---------------------------------------% */

	if (*rvec) {
	    zscal_(ncv, &rnorm, &workl[ihbds], &c__1);
	}

	i__1 = *ncv;
	for (k = 1; k <= i__1; ++k) {
	    i__2 = iheig + k - 1;
	    temp.r = workl[i__2].r, temp.i = workl[i__2].i;
	    i__2 = ihbds + k - 1;
	    z_div(&z__2, &workl[ihbds + k - 1], &temp);
	    z_div(&z__1, &z__2, &temp);
	    workl[i__2].r = z__1.r, workl[i__2].i = z__1.i;
/* L50: */
	}

    }

/*     %-----------------------------------------------------------% */
/*     | *  Transform the Ritz values back to the original system. | */
/*     |    For TYPE = 'SHIFTI' the transformation is              | */
/*     |             lambda = 1/theta + sigma                      | */
/*     | NOTES:                                                    | */
/*     | *The Ritz vectors are not affected by the transformation. | */
/*     %-----------------------------------------------------------% */

    if (s_cmp(type__, "SHIFTI", (ftnlen)6, (ftnlen)6) == 0) {
	i__1 = nconv;
	for (k = 1; k <= i__1; ++k) {
	    i__2 = k;
	    z_div(&z__2, &c_b1, &workl[iheig + k - 1]);
	    z__1.r = z__2.r + sigma->r, z__1.i = z__2.i + sigma->i;
	    d__[i__2].r = z__1.r, d__[i__2].i = z__1.i;
/* L60: */
	}
    }

    if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) != 0 && msglvl > 1) {
	zvout_(&debug_1.logfil, &nconv, &d__[1], &debug_1.ndigit, "_neupd: U"
		"ntransformed Ritz values.", (ftnlen)34);
	zvout_(&debug_1.logfil, &nconv, &workl[ihbds], &debug_1.ndigit, "_ne"
		"upd: Ritz estimates of the untransformed Ritz values.", (
		ftnlen)56);
    } else if (msglvl > 1) {
	zvout_(&debug_1.logfil, &nconv, &d__[1], &debug_1.ndigit, "_neupd: C"
		"onverged Ritz values.", (ftnlen)30);
	zvout_(&debug_1.logfil, &nconv, &workl[ihbds], &debug_1.ndigit, "_ne"
		"upd: Associated Ritz estimates.", (ftnlen)34);
    }

/*     %-------------------------------------------------% */
/*     | Eigenvector Purification step. Formally perform | */
/*     | one of inverse subspace iteration. Only used    | */
/*     | for MODE = 3. See reference 3.                  | */
/*     %-------------------------------------------------% */

    if (*rvec && *(unsigned char *)howmny == 'A' && s_cmp(type__, "SHIFTI", (
	    ftnlen)6, (ftnlen)6) == 0) {

/*        %------------------------------------------------% */
/*        | Purify the computed Ritz vectors by adding a   | */
/*        | little bit of the residual vector:             | */
/*        |                      T                         | */
/*        |          resid(:)*( e    s ) / theta           | */
/*        |                      NCV                       | */
/*        | where H s = s theta.                           | */
/*        %------------------------------------------------% */

	i__1 = nconv;
	for (j = 1; j <= i__1; ++j) {
	    i__2 = iheig + j - 1;
	    if (workl[i__2].r != 0. || workl[i__2].i != 0.) {
		i__2 = j;
		z_div(&z__1, &workl[invsub + (j - 1) * ldq + *ncv - 1], &
			workl[iheig + j - 1]);
		workev[i__2].r = z__1.r, workev[i__2].i = z__1.i;
	    }
/* L100: */
	}
/*        %---------------------------------------% */
/*        | Perform a rank one update to Z and    | */
/*        | purify all the Ritz vectors together. | */
/*        %---------------------------------------% */

	zgeru_(n, &nconv, &c_b1, &resid[1], &c__1, &workev[1], &c__1, &z__[
		z_offset], ldz);

    }

L9000:

    return 0;

/*     %---------------% */
/*     | End of zneupd | */
/*     %---------------% */

} /* zneupd_ */