/*  -- 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 integer c__1 = 1;
static logical c_true = TRUE_;
static logical c_false = FALSE_;

/* > \brief \b DTRSNA   

    =========== DOCUMENTATION ===========   

   Online html documentation available at   
              http://www.netlib.org/lapack/explore-html/   

   > \htmlonly   
   > Download DTRSNA + dependencies   
   > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna.
f">   
   > [TGZ]</a>   
   > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna.
f">   
   > [ZIP]</a>   
   > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna.
f">   
   > [TXT]</a>   
   > \endhtmlonly   

    Definition:   
    ===========   

         SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,   
                            LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,   
                            INFO )   

         CHARACTER          HOWMNY, JOB   
         INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N   
         LOGICAL            SELECT( * )   
         INTEGER            IWORK( * )   
         DOUBLE PRECISION   S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),   
        $                   VR( LDVR, * ), WORK( LDWORK, * )   


   > \par Purpose:   
    =============   
   >   
   > \verbatim   
   >   
   > DTRSNA estimates reciprocal condition numbers for specified   
   > eigenvalues and/or right eigenvectors of a real upper   
   > quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q   
   > orthogonal).   
   >   
   > T must be in Schur canonical form (as returned by DHSEQR), that is,   
   > block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each   
   > 2-by-2 diagonal block has its diagonal elements equal and its   
   > off-diagonal elements of opposite sign.   
   > \endverbatim   

    Arguments:   
    ==========   

   > \param[in] JOB   
   > \verbatim   
   >          JOB is CHARACTER*1   
   >          Specifies whether condition numbers are required for   
   >          eigenvalues (S) or eigenvectors (SEP):   
   >          = 'E': for eigenvalues only (S);   
   >          = 'V': for eigenvectors only (SEP);   
   >          = 'B': for both eigenvalues and eigenvectors (S and SEP).   
   > \endverbatim   
   >   
   > \param[in] HOWMNY   
   > \verbatim   
   >          HOWMNY is CHARACTER*1   
   >          = 'A': compute condition numbers for all eigenpairs;   
   >          = 'S': compute condition numbers for selected eigenpairs   
   >                 specified by the array SELECT.   
   > \endverbatim   
   >   
   > \param[in] SELECT   
   > \verbatim   
   >          SELECT is LOGICAL array, dimension (N)   
   >          If HOWMNY = 'S', SELECT specifies the eigenpairs for which   
   >          condition numbers are required. To select condition numbers   
   >          for the eigenpair corresponding to a real eigenvalue w(j),   
   >          SELECT(j) must be set to .TRUE.. To select condition numbers   
   >          corresponding to a complex conjugate pair of eigenvalues w(j)   
   >          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be   
   >          set to .TRUE..   
   >          If HOWMNY = 'A', SELECT is not referenced.   
   > \endverbatim   
   >   
   > \param[in] N   
   > \verbatim   
   >          N is INTEGER   
   >          The order of the matrix T. N >= 0.   
   > \endverbatim   
   >   
   > \param[in] T   
   > \verbatim   
   >          T is DOUBLE PRECISION array, dimension (LDT,N)   
   >          The upper quasi-triangular matrix T, in Schur canonical form.   
   > \endverbatim   
   >   
   > \param[in] LDT   
   > \verbatim   
   >          LDT is INTEGER   
   >          The leading dimension of the array T. LDT >= max(1,N).   
   > \endverbatim   
   >   
   > \param[in] VL   
   > \verbatim   
   >          VL is DOUBLE PRECISION array, dimension (LDVL,M)   
   >          If JOB = 'E' or 'B', VL must contain left eigenvectors of T   
   >          (or of any Q*T*Q**T with Q orthogonal), corresponding to the   
   >          eigenpairs specified by HOWMNY and SELECT. The eigenvectors   
   >          must be stored in consecutive columns of VL, as returned by   
   >          DHSEIN or DTREVC.   
   >          If JOB = 'V', VL is not referenced.   
   > \endverbatim   
   >   
   > \param[in] LDVL   
   > \verbatim   
   >          LDVL is INTEGER   
   >          The leading dimension of the array VL.   
   >          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.   
   > \endverbatim   
   >   
   > \param[in] VR   
   > \verbatim   
   >          VR is DOUBLE PRECISION array, dimension (LDVR,M)   
   >          If JOB = 'E' or 'B', VR must contain right eigenvectors of T   
   >          (or of any Q*T*Q**T with Q orthogonal), corresponding to the   
   >          eigenpairs specified by HOWMNY and SELECT. The eigenvectors   
   >          must be stored in consecutive columns of VR, as returned by   
   >          DHSEIN or DTREVC.   
   >          If JOB = 'V', VR is not referenced.   
   > \endverbatim   
   >   
   > \param[in] LDVR   
   > \verbatim   
   >          LDVR is INTEGER   
   >          The leading dimension of the array VR.   
   >          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.   
   > \endverbatim   
   >   
   > \param[out] S   
   > \verbatim   
   >          S is DOUBLE PRECISION array, dimension (MM)   
   >          If JOB = 'E' or 'B', the reciprocal condition numbers of the   
   >          selected eigenvalues, stored in consecutive elements of the   
   >          array. For a complex conjugate pair of eigenvalues two   
   >          consecutive elements of S are set to the same value. Thus   
   >          S(j), SEP(j), and the j-th columns of VL and VR all   
   >          correspond to the same eigenpair (but not in general the   
   >          j-th eigenpair, unless all eigenpairs are selected).   
   >          If JOB = 'V', S is not referenced.   
   > \endverbatim   
   >   
   > \param[out] SEP   
   > \verbatim   
   >          SEP is DOUBLE PRECISION array, dimension (MM)   
   >          If JOB = 'V' or 'B', the estimated reciprocal condition   
   >          numbers of the selected eigenvectors, stored in consecutive   
   >          elements of the array. For a complex eigenvector two   
   >          consecutive elements of SEP are set to the same value. If   
   >          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)   
   >          is set to 0; this can only occur when the true value would be   
   >          very small anyway.   
   >          If JOB = 'E', SEP is not referenced.   
   > \endverbatim   
   >   
   > \param[in] MM   
   > \verbatim   
   >          MM is INTEGER   
   >          The number of elements in the arrays S (if JOB = 'E' or 'B')   
   >           and/or SEP (if JOB = 'V' or 'B'). MM >= M.   
   > \endverbatim   
   >   
   > \param[out] M   
   > \verbatim   
   >          M is INTEGER   
   >          The number of elements of the arrays S and/or SEP actually   
   >          used to store the estimated condition numbers.   
   >          If HOWMNY = 'A', M is set to N.   
   > \endverbatim   
   >   
   > \param[out] WORK   
   > \verbatim   
   >          WORK is DOUBLE PRECISION array, dimension (LDWORK,N+6)   
   >          If JOB = 'E', WORK is not referenced.   
   > \endverbatim   
   >   
   > \param[in] LDWORK   
   > \verbatim   
   >          LDWORK is INTEGER   
   >          The leading dimension of the array WORK.   
   >          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.   
   > \endverbatim   
   >   
   > \param[out] IWORK   
   > \verbatim   
   >          IWORK is INTEGER array, dimension (2*(N-1))   
   >          If JOB = 'E', IWORK is not referenced.   
   > \endverbatim   
   >   
   > \param[out] INFO   
   > \verbatim   
   >          INFO is INTEGER   
   >          = 0: successful exit   
   >          < 0: if INFO = -i, the i-th argument had an illegal value   
   > \endverbatim   

    Authors:   
    ========   

   > \author Univ. of Tennessee   
   > \author Univ. of California Berkeley   
   > \author Univ. of Colorado Denver   
   > \author NAG Ltd.   

   > \date November 2011   

   > \ingroup doubleOTHERcomputational   

   > \par Further Details:   
    =====================   
   >   
   > \verbatim   
   >   
   >  The reciprocal of the condition number of an eigenvalue lambda is   
   >  defined as   
   >   
   >          S(lambda) = |v**T*u| / (norm(u)*norm(v))   
   >   
   >  where u and v are the right and left eigenvectors of T corresponding   
   >  to lambda; v**T denotes the transpose of v, and norm(u)   
   >  denotes the Euclidean norm. These reciprocal condition numbers always   
   >  lie between zero (very badly conditioned) and one (very well   
   >  conditioned). If n = 1, S(lambda) is defined to be 1.   
   >   
   >  An approximate error bound for a computed eigenvalue W(i) is given by   
   >   
   >                      EPS * norm(T) / S(i)   
   >   
   >  where EPS is the machine precision.   
   >   
   >  The reciprocal of the condition number of the right eigenvector u   
   >  corresponding to lambda is defined as follows. Suppose   
   >   
   >              T = ( lambda  c  )   
   >                  (   0    T22 )   
   >   
   >  Then the reciprocal condition number is   
   >   
   >          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )   
   >   
   >  where sigma-min denotes the smallest singular value. We approximate   
   >  the smallest singular value by the reciprocal of an estimate of the   
   >  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is   
   >  defined to be abs(T(1,1)).   
   >   
   >  An approximate error bound for a computed right eigenvector VR(i)   
   >  is given by   
   >   
   >                      EPS * norm(T) / SEP(i)   
   > \endverbatim   
   >   
    =====================================================================   
   Subroutine */ int igraphdtrsna_(char *job, char *howmny, logical *select, 
	integer *n, doublereal *t, integer *ldt, doublereal *vl, integer *
	ldvl, doublereal *vr, integer *ldvr, doublereal *s, doublereal *sep, 
	integer *mm, integer *m, doublereal *work, integer *ldwork, integer *
	iwork, integer *info)
{
    /* System generated locals */
    integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, 
	    work_dim1, work_offset, i__1, i__2;
    doublereal d__1, d__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    integer i__, j, k, n2;
    doublereal cs;
    integer nn, ks;
    doublereal sn, mu, eps, est;
    integer kase;
    doublereal cond;
    extern doublereal igraphddot_(integer *, doublereal *, integer *, doublereal *, 
	    integer *);
    logical pair;
    integer ierr;
    doublereal dumm, prod;
    integer ifst;
    doublereal lnrm;
    integer ilst;
    doublereal rnrm;
    extern doublereal igraphdnrm2_(integer *, doublereal *, integer *);
    doublereal prod1, prod2, scale, delta;
    extern logical igraphlsame_(char *, char *);
    integer isave[3];
    logical wants;
    doublereal dummy[1];
    extern /* Subroutine */ int igraphdlacn2_(integer *, doublereal *, doublereal *,
	     integer *, doublereal *, integer *, integer *);
    extern doublereal igraphdlapy2_(doublereal *, doublereal *);
    extern /* Subroutine */ int igraphdlabad_(doublereal *, doublereal *);
    extern doublereal igraphdlamch_(char *);
    extern /* Subroutine */ int igraphdlacpy_(char *, integer *, integer *, 
	    doublereal *, integer *, doublereal *, integer *), 
	    igraphxerbla_(char *, integer *, ftnlen);
    doublereal bignum;
    logical wantbh;
    extern /* Subroutine */ int igraphdlaqtr_(logical *, logical *, integer *, 
	    doublereal *, integer *, doublereal *, doublereal *, doublereal *,
	     doublereal *, doublereal *, integer *), igraphdtrexc_(char *, integer *
	    , doublereal *, integer *, doublereal *, integer *, integer *, 
	    integer *, doublereal *, integer *);
    logical somcon;
    doublereal smlnum;
    logical wantsp;


/*  -- LAPACK computational routine (version 3.4.0) --   
    -- LAPACK is a software package provided by Univ. of Tennessee,    --   
    -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--   
       November 2011   


    =====================================================================   


       Decode and test the input parameters   

       Parameter adjustments */
    --select;
    t_dim1 = *ldt;
    t_offset = 1 + t_dim1;
    t -= t_offset;
    vl_dim1 = *ldvl;
    vl_offset = 1 + vl_dim1;
    vl -= vl_offset;
    vr_dim1 = *ldvr;
    vr_offset = 1 + vr_dim1;
    vr -= vr_offset;
    --s;
    --sep;
    work_dim1 = *ldwork;
    work_offset = 1 + work_dim1;
    work -= work_offset;
    --iwork;

    /* Function Body */
    wantbh = igraphlsame_(job, "B");
    wants = igraphlsame_(job, "E") || wantbh;
    wantsp = igraphlsame_(job, "V") || wantbh;

    somcon = igraphlsame_(howmny, "S");

    *info = 0;
    if (! wants && ! wantsp) {
	*info = -1;
    } else if (! igraphlsame_(howmny, "A") && ! somcon) {
	*info = -2;
    } else if (*n < 0) {
	*info = -4;
    } else if (*ldt < max(1,*n)) {
	*info = -6;
    } else if (*ldvl < 1 || wants && *ldvl < *n) {
	*info = -8;
    } else if (*ldvr < 1 || wants && *ldvr < *n) {
	*info = -10;
    } else {

/*        Set M to the number of eigenpairs for which condition numbers   
          are required, and test MM. */

	if (somcon) {
	    *m = 0;
	    pair = FALSE_;
	    i__1 = *n;
	    for (k = 1; k <= i__1; ++k) {
		if (pair) {
		    pair = FALSE_;
		} else {
		    if (k < *n) {
			if (t[k + 1 + k * t_dim1] == 0.) {
			    if (select[k]) {
				++(*m);
			    }
			} else {
			    pair = TRUE_;
			    if (select[k] || select[k + 1]) {
				*m += 2;
			    }
			}
		    } else {
			if (select[*n]) {
			    ++(*m);
			}
		    }
		}
/* L10: */
	    }
	} else {
	    *m = *n;
	}

	if (*mm < *m) {
	    *info = -13;
	} else if (*ldwork < 1 || wantsp && *ldwork < *n) {
	    *info = -16;
	}
    }
    if (*info != 0) {
	i__1 = -(*info);
	igraphxerbla_("DTRSNA", &i__1, (ftnlen)6);
	return 0;
    }

/*     Quick return if possible */

    if (*n == 0) {
	return 0;
    }

    if (*n == 1) {
	if (somcon) {
	    if (! select[1]) {
		return 0;
	    }
	}
	if (wants) {
	    s[1] = 1.;
	}
	if (wantsp) {
	    sep[1] = (d__1 = t[t_dim1 + 1], abs(d__1));
	}
	return 0;
    }

/*     Get machine constants */

    eps = igraphdlamch_("P");
    smlnum = igraphdlamch_("S") / eps;
    bignum = 1. / smlnum;
    igraphdlabad_(&smlnum, &bignum);

    ks = 0;
    pair = FALSE_;
    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {

/*        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block. */

	if (pair) {
	    pair = FALSE_;
	    goto L60;
	} else {
	    if (k < *n) {
		pair = t[k + 1 + k * t_dim1] != 0.;
	    }
	}

/*        Determine whether condition numbers are required for the k-th   
          eigenpair. */

	if (somcon) {
	    if (pair) {
		if (! select[k] && ! select[k + 1]) {
		    goto L60;
		}
	    } else {
		if (! select[k]) {
		    goto L60;
		}
	    }
	}

	++ks;

	if (wants) {

/*           Compute the reciprocal condition number of the k-th   
             eigenvalue. */

	    if (! pair) {

/*              Real eigenvalue. */

		prod = igraphddot_(n, &vr[ks * vr_dim1 + 1], &c__1, &vl[ks * 
			vl_dim1 + 1], &c__1);
		rnrm = igraphdnrm2_(n, &vr[ks * vr_dim1 + 1], &c__1);
		lnrm = igraphdnrm2_(n, &vl[ks * vl_dim1 + 1], &c__1);
		s[ks] = abs(prod) / (rnrm * lnrm);
	    } else {

/*              Complex eigenvalue. */

		prod1 = igraphddot_(n, &vr[ks * vr_dim1 + 1], &c__1, &vl[ks * 
			vl_dim1 + 1], &c__1);
		prod1 += igraphddot_(n, &vr[(ks + 1) * vr_dim1 + 1], &c__1, &vl[(ks 
			+ 1) * vl_dim1 + 1], &c__1);
		prod2 = igraphddot_(n, &vl[ks * vl_dim1 + 1], &c__1, &vr[(ks + 1) * 
			vr_dim1 + 1], &c__1);
		prod2 -= igraphddot_(n, &vl[(ks + 1) * vl_dim1 + 1], &c__1, &vr[ks *
			 vr_dim1 + 1], &c__1);
		d__1 = igraphdnrm2_(n, &vr[ks * vr_dim1 + 1], &c__1);
		d__2 = igraphdnrm2_(n, &vr[(ks + 1) * vr_dim1 + 1], &c__1);
		rnrm = igraphdlapy2_(&d__1, &d__2);
		d__1 = igraphdnrm2_(n, &vl[ks * vl_dim1 + 1], &c__1);
		d__2 = igraphdnrm2_(n, &vl[(ks + 1) * vl_dim1 + 1], &c__1);
		lnrm = igraphdlapy2_(&d__1, &d__2);
		cond = igraphdlapy2_(&prod1, &prod2) / (rnrm * lnrm);
		s[ks] = cond;
		s[ks + 1] = cond;
	    }
	}

	if (wantsp) {

/*           Estimate the reciprocal condition number of the k-th   
             eigenvector.   

             Copy the matrix T to the array WORK and swap the diagonal   
             block beginning at T(k,k) to the (1,1) position. */

	    igraphdlacpy_("Full", n, n, &t[t_offset], ldt, &work[work_offset], 
		    ldwork);
	    ifst = k;
	    ilst = 1;
	    igraphdtrexc_("No Q", n, &work[work_offset], ldwork, dummy, &c__1, &
		    ifst, &ilst, &work[(*n + 1) * work_dim1 + 1], &ierr);

	    if (ierr == 1 || ierr == 2) {

/*              Could not swap because blocks not well separated */

		scale = 1.;
		est = bignum;
	    } else {

/*              Reordering successful */

		if (work[work_dim1 + 2] == 0.) {

/*                 Form C = T22 - lambda*I in WORK(2:N,2:N). */

		    i__2 = *n;
		    for (i__ = 2; i__ <= i__2; ++i__) {
			work[i__ + i__ * work_dim1] -= work[work_dim1 + 1];
/* L20: */
		    }
		    n2 = 1;
		    nn = *n - 1;
		} else {

/*                 Triangularize the 2 by 2 block by unitary   
                   transformation U = [  cs   i*ss ]   
                                      [ i*ss   cs  ].   
                   such that the (1,1) position of WORK is complex   
                   eigenvalue lambda with positive imaginary part. (2,2)   
                   position of WORK is the complex eigenvalue lambda   
                   with negative imaginary  part. */

		    mu = sqrt((d__1 = work[(work_dim1 << 1) + 1], abs(d__1))) 
			    * sqrt((d__2 = work[work_dim1 + 2], abs(d__2)));
		    delta = igraphdlapy2_(&mu, &work[work_dim1 + 2]);
		    cs = mu / delta;
		    sn = -work[work_dim1 + 2] / delta;

/*                 Form   

                   C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]   
                                            [   mu                     ]   
                                            [         ..               ]   
                                            [             ..           ]   
                                            [                  mu      ]   
                   where C**T is transpose of matrix C,   
                   and RWORK is stored starting in the N+1-st column of   
                   WORK. */

		    i__2 = *n;
		    for (j = 3; j <= i__2; ++j) {
			work[j * work_dim1 + 2] = cs * work[j * work_dim1 + 2]
				;
			work[j + j * work_dim1] -= work[work_dim1 + 1];
/* L30: */
		    }
		    work[(work_dim1 << 1) + 2] = 0.;

		    work[(*n + 1) * work_dim1 + 1] = mu * 2.;
		    i__2 = *n - 1;
		    for (i__ = 2; i__ <= i__2; ++i__) {
			work[i__ + (*n + 1) * work_dim1] = sn * work[(i__ + 1)
				 * work_dim1 + 1];
/* L40: */
		    }
		    n2 = 2;
		    nn = *n - 1 << 1;
		}

/*              Estimate norm(inv(C**T)) */

		est = 0.;
		kase = 0;
L50:
		igraphdlacn2_(&nn, &work[(*n + 2) * work_dim1 + 1], &work[(*n + 4) *
			 work_dim1 + 1], &iwork[1], &est, &kase, isave);
		if (kase != 0) {
		    if (kase == 1) {
			if (n2 == 1) {

/*                       Real eigenvalue: solve C**T*x = scale*c. */

			    i__2 = *n - 1;
			    igraphdlaqtr_(&c_true, &c_true, &i__2, &work[(work_dim1 
				    << 1) + 2], ldwork, dummy, &dumm, &scale, 
				    &work[(*n + 4) * work_dim1 + 1], &work[(*
				    n + 6) * work_dim1 + 1], &ierr);
			} else {

/*                       Complex eigenvalue: solve   
                         C**T*(p+iq) = scale*(c+id) in real arithmetic. */

			    i__2 = *n - 1;
			    igraphdlaqtr_(&c_true, &c_false, &i__2, &work[(
				    work_dim1 << 1) + 2], ldwork, &work[(*n + 
				    1) * work_dim1 + 1], &mu, &scale, &work[(*
				    n + 4) * work_dim1 + 1], &work[(*n + 6) * 
				    work_dim1 + 1], &ierr);
			}
		    } else {
			if (n2 == 1) {

/*                       Real eigenvalue: solve C*x = scale*c. */

			    i__2 = *n - 1;
			    igraphdlaqtr_(&c_false, &c_true, &i__2, &work[(
				    work_dim1 << 1) + 2], ldwork, dummy, &
				    dumm, &scale, &work[(*n + 4) * work_dim1 
				    + 1], &work[(*n + 6) * work_dim1 + 1], &
				    ierr);
			} else {

/*                       Complex eigenvalue: solve   
                         C*(p+iq) = scale*(c+id) in real arithmetic. */

			    i__2 = *n - 1;
			    igraphdlaqtr_(&c_false, &c_false, &i__2, &work[(
				    work_dim1 << 1) + 2], ldwork, &work[(*n + 
				    1) * work_dim1 + 1], &mu, &scale, &work[(*
				    n + 4) * work_dim1 + 1], &work[(*n + 6) * 
				    work_dim1 + 1], &ierr);

			}
		    }

		    goto L50;
		}
	    }

	    sep[ks] = scale / max(est,smlnum);
	    if (pair) {
		sep[ks + 1] = sep[ks];
	    }
	}

	if (pair) {
	    ++ks;
	}

L60:
	;
    }
    return 0;

/*     End of DTRSNA */

} /* igraphdtrsna_ */