/*  -- translated by f2c (version 20100827).
   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"

/* Table of constant values */

static logical c_true = TRUE_;
static integer c__1 = 1;
static doublereal c_b18 = 1.;
static doublereal c_b20 = 0.;

/* -----------------------------------------------------------------------   
   \BeginDoc   

   \Name: dneigh   

   \Description:   
    Compute the eigenvalues of the current upper Hessenberg matrix   
    and the corresponding Ritz estimates given the current residual norm.   

   \Usage:   
    call dneigh   
       ( RNORM, N, H, LDH, RITZR, RITZI, BOUNDS, Q, LDQ, WORKL, IERR )   

   \Arguments   
    RNORM   Double precision scalar.  (INPUT)   
            Residual norm corresponding to the current upper Hessenberg   
            matrix H.   

    N       Integer.  (INPUT)   
            Size of the matrix H.   

    H       Double precision N by N array.  (INPUT)   
            H contains the current upper Hessenberg matrix.   

    LDH     Integer.  (INPUT)   
            Leading dimension of H exactly as declared in the calling   
            program.   

    RITZR,  Double precision arrays of length N.  (OUTPUT)   
    RITZI   On output, RITZR(1:N) (resp. RITZI(1:N)) contains the real   
            (respectively imaginary) parts of the eigenvalues of H.   

    BOUNDS  Double precision array of length N.  (OUTPUT)   
            On output, BOUNDS contains the Ritz estimates associated with   
            the eigenvalues RITZR and RITZI.  This is equal to RNORM   
            times the last components of the eigenvectors corresponding   
            to the eigenvalues in RITZR and RITZI.   

    Q       Double precision N by N array.  (WORKSPACE)   
            Workspace needed to store the eigenvectors of H.   

    LDQ     Integer.  (INPUT)   
            Leading dimension of Q exactly as declared in the calling   
            program.   

    WORKL   Double precision work array of length N**2 + 3*N.  (WORKSPACE)   
            Private (replicated) array on each PE or array allocated on   
            the front end.  This is needed to keep the full Schur form   
            of H and also in the calculation of the eigenvectors of H.   

    IERR    Integer.  (OUTPUT)   
            Error exit flag from dlaqrb or dtrevc.   

   \EndDoc   

   -----------------------------------------------------------------------   

   \BeginLib   

   \Local variables:   
       xxxxxx  real   

   \Routines called:   
       dlaqrb  ARPACK routine to compute the real Schur form of an   
               upper Hessenberg matrix and last row of the Schur vectors.   
       second  ARPACK utility routine for timing.   
       dmout   ARPACK utility routine that prints matrices   
       dvout   ARPACK utility routine that prints vectors.   
       dlacpy  LAPACK matrix copy routine.   
       dlapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.   
       dtrevc  LAPACK routine to compute the eigenvectors of a matrix   
               in upper quasi-triangular form   
       dgemv   Level 2 BLAS routine for matrix vector multiplication.   
       dcopy   Level 1 BLAS that copies one vector to another .   
       dnrm2   Level 1 BLAS that computes the norm of a vector.   
       dscal   Level 1 BLAS that scales a vector.   


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

   \Revision history:   
       xx/xx/92: Version ' 2.1'   

   \SCCS Information: @(#)   
   FILE: neigh.F   SID: 2.3   DATE OF SID: 4/20/96   RELEASE: 2   

   \Remarks   
       None   

   \EndLib   

   -----------------------------------------------------------------------   

   Subroutine */ int igraphdneigh_(doublereal *rnorm, integer *n, doublereal *h__, 
	integer *ldh, doublereal *ritzr, doublereal *ritzi, doublereal *
	bounds, doublereal *q, integer *ldq, doublereal *workl, integer *ierr)
{
    /* System generated locals */
    integer h_dim1, h_offset, q_dim1, q_offset, i__1;
    doublereal d__1, d__2;

    /* Local variables */
    integer i__;
    real t0, t1;
    doublereal vl[1], temp;
    extern doublereal igraphdnrm2_(integer *, doublereal *, integer *);
    extern /* Subroutine */ int igraphdscal_(integer *, doublereal *, doublereal *, 
	    integer *);
    integer iconj;
    extern /* Subroutine */ int igraphdgemv_(char *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, doublereal *, integer *), igraphdmout_(integer *, 
	    integer *, integer *, doublereal *, integer *, integer *, char *, 
	    ftnlen), igraphdvout_(integer *, integer *, doublereal *, integer *, 
	    char *, ftnlen);
    extern doublereal igraphdlapy2_(doublereal *, doublereal *);
    extern /* Subroutine */ int igraphdlaqrb_(logical *, integer *, integer *, 
	    integer *, doublereal *, integer *, doublereal *, doublereal *, 
	    doublereal *, integer *);
    integer mneigh = 0;
    extern /* Subroutine */ int igraphsecond_(real *), igraphdlacpy_(char *, integer *, 
	    integer *, doublereal *, integer *, doublereal *, integer *);
    integer logfil, ndigit;
    logical select[1];
    real tneigh = 0.;
    extern /* Subroutine */ int igraphdtrevc_(char *, char *, logical *, integer *, 
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
	    integer *, integer *, integer *, doublereal *, integer *);
    integer msglvl;


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


       %------------------%   
       | Scalar Arguments |   
       %------------------%   


       %-----------------%   
       | Array Arguments |   
       %-----------------%   


       %------------%   
       | Parameters |   
       %------------%   


       %------------------------%   
       | Local Scalars & Arrays |   
       %------------------------%   


       %----------------------%   
       | External Subroutines |   
       %----------------------%   


       %--------------------%   
       | External Functions |   
       %--------------------%   


       %---------------------%   
       | Intrinsic Functions |   
       %---------------------%   


       %-----------------------%   
       | Executable Statements |   
       %-----------------------%   


       %-------------------------------%   
       | Initialize timing statistics  |   
       | & message level for debugging |   
       %-------------------------------%   

       Parameter adjustments */
    --workl;
    --bounds;
    --ritzi;
    --ritzr;
    h_dim1 = *ldh;
    h_offset = 1 + h_dim1;
    h__ -= h_offset;
    q_dim1 = *ldq;
    q_offset = 1 + q_dim1;
    q -= q_offset;

    /* Function Body */
    igraphsecond_(&t0);
    msglvl = mneigh;

    if (msglvl > 2) {
	igraphdmout_(&logfil, n, n, &h__[h_offset], ldh, &ndigit, "_neigh: Enterin"
		"g upper Hessenberg matrix H ", (ftnlen)43);
    }

/*     %-----------------------------------------------------------%   
       | 1. Compute the eigenvalues, the last components of the    |   
       |    corresponding Schur vectors and the full Schur form T  |   
       |    of the current upper Hessenberg matrix H.              |   
       | dlaqrb returns the full Schur form of H in WORKL(1:N**2)  |   
       | and the last components of the Schur vectors in BOUNDS.   |   
       %-----------------------------------------------------------% */

    igraphdlacpy_("All", n, n, &h__[h_offset], ldh, &workl[1], n);
    igraphdlaqrb_(&c_true, n, &c__1, n, &workl[1], n, &ritzr[1], &ritzi[1], &bounds[
	    1], ierr);
    if (*ierr != 0) {
	goto L9000;
    }

    if (msglvl > 1) {
	igraphdvout_(&logfil, n, &bounds[1], &ndigit, "_neigh: last row of the Sch"
		"ur matrix for H", (ftnlen)42);
    }

/*     %-----------------------------------------------------------%   
       | 2. Compute the eigenvectors of the full Schur form T and  |   
       |    apply the last components of the Schur vectors to get  |   
       |    the last components of the corresponding eigenvectors. |   
       | Remember that if the i-th and (i+1)-st eigenvalues are    |   
       | complex conjugate pairs, then the real & imaginary part   |   
       | of the eigenvector components are split across adjacent   |   
       | columns of Q.                                             |   
       %-----------------------------------------------------------% */

    igraphdtrevc_("R", "A", select, n, &workl[1], n, vl, n, &q[q_offset], ldq, n, n,
	     &workl[*n * *n + 1], ierr);

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

/*     %------------------------------------------------%   
       | Scale the returning eigenvectors so that their |   
       | euclidean norms are all one. LAPACK subroutine |   
       | dtrevc returns each eigenvector normalized so  |   
       | that the element of largest magnitude has      |   
       | magnitude 1; here the magnitude of a complex   |   
       | number (x,y) is taken to be |x| + |y|.         |   
       %------------------------------------------------% */

    iconj = 0;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	if ((d__1 = ritzi[i__], abs(d__1)) <= 0.) {

/*           %----------------------%   
             | Real eigenvalue case |   
             %----------------------% */

	    temp = igraphdnrm2_(n, &q[i__ * q_dim1 + 1], &c__1);
	    d__1 = 1. / temp;
	    igraphdscal_(n, &d__1, &q[i__ * q_dim1 + 1], &c__1);
	} else {

/*           %-------------------------------------------%   
             | Complex conjugate pair case. Note that    |   
             | since the real and imaginary part of      |   
             | the eigenvector are stored in consecutive |   
             | columns, we further normalize by the      |   
             | square root of two.                       |   
             %-------------------------------------------% */

	    if (iconj == 0) {
		d__1 = igraphdnrm2_(n, &q[i__ * q_dim1 + 1], &c__1);
		d__2 = igraphdnrm2_(n, &q[(i__ + 1) * q_dim1 + 1], &c__1);
		temp = igraphdlapy2_(&d__1, &d__2);
		d__1 = 1. / temp;
		igraphdscal_(n, &d__1, &q[i__ * q_dim1 + 1], &c__1);
		d__1 = 1. / temp;
		igraphdscal_(n, &d__1, &q[(i__ + 1) * q_dim1 + 1], &c__1);
		iconj = 1;
	    } else {
		iconj = 0;
	    }
	}
/* L10: */
    }

    igraphdgemv_("T", n, n, &c_b18, &q[q_offset], ldq, &bounds[1], &c__1, &c_b20, &
	    workl[1], &c__1);

    if (msglvl > 1) {
	igraphdvout_(&logfil, n, &workl[1], &ndigit, "_neigh: Last row of the eige"
		"nvector matrix for H", (ftnlen)48);
    }

/*     %----------------------------%   
       | Compute the Ritz estimates |   
       %----------------------------% */

    iconj = 0;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	if ((d__1 = ritzi[i__], abs(d__1)) <= 0.) {

/*           %----------------------%   
             | Real eigenvalue case |   
             %----------------------% */

	    bounds[i__] = *rnorm * (d__1 = workl[i__], abs(d__1));
	} else {

/*           %-------------------------------------------%   
             | Complex conjugate pair case. Note that    |   
             | since the real and imaginary part of      |   
             | the eigenvector are stored in consecutive |   
             | columns, we need to take the magnitude    |   
             | of the last components of the two vectors |   
             %-------------------------------------------% */

	    if (iconj == 0) {
		bounds[i__] = *rnorm * igraphdlapy2_(&workl[i__], &workl[i__ + 1]);
		bounds[i__ + 1] = bounds[i__];
		iconj = 1;
	    } else {
		iconj = 0;
	    }
	}
/* L20: */
    }

    if (msglvl > 2) {
	igraphdvout_(&logfil, n, &ritzr[1], &ndigit, "_neigh: Real part of the eig"
		"envalues of H", (ftnlen)41);
	igraphdvout_(&logfil, n, &ritzi[1], &ndigit, "_neigh: Imaginary part of th"
		"e eigenvalues of H", (ftnlen)46);
	igraphdvout_(&logfil, n, &bounds[1], &ndigit, "_neigh: Ritz estimates for "
		"the eigenvalues of H", (ftnlen)47);
    }

    igraphsecond_(&t1);
    tneigh += t1 - t0;

L9000:
    return 0;

/*     %---------------%   
       | End of dneigh |   
       %---------------% */

} /* igraphdneigh_ */