/*  -- 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_n1 = -1;

/* > \brief \b DTRSEN   

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

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

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

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

         SUBROUTINE DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,   
                            M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )   

         CHARACTER          COMPQ, JOB   
         INTEGER            INFO, LDQ, LDT, LIWORK, LWORK, M, N   
         DOUBLE PRECISION   S, SEP   
         LOGICAL            SELECT( * )   
         INTEGER            IWORK( * )   
         DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),   
        $                   WR( * )   


   > \par Purpose:   
    =============   
   >   
   > \verbatim   
   >   
   > DTRSEN reorders the real Schur factorization of a real matrix   
   > A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in   
   > the leading diagonal blocks of the upper quasi-triangular matrix T,   
   > and the leading columns of Q form an orthonormal basis of the   
   > corresponding right invariant subspace.   
   >   
   > Optionally the routine computes the reciprocal condition numbers of   
   > the cluster of eigenvalues and/or the invariant subspace.   
   >   
   > 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 the   
   >          cluster of eigenvalues (S) or the invariant subspace (SEP):   
   >          = 'N': none;   
   >          = 'E': for eigenvalues only (S);   
   >          = 'V': for invariant subspace only (SEP);   
   >          = 'B': for both eigenvalues and invariant subspace (S and   
   >                 SEP).   
   > \endverbatim   
   >   
   > \param[in] COMPQ   
   > \verbatim   
   >          COMPQ is CHARACTER*1   
   >          = 'V': update the matrix Q of Schur vectors;   
   >          = 'N': do not update Q.   
   > \endverbatim   
   >   
   > \param[in] SELECT   
   > \verbatim   
   >          SELECT is LOGICAL array, dimension (N)   
   >          SELECT specifies the eigenvalues in the selected cluster. To   
   >          select a real eigenvalue w(j), SELECT(j) must be set to   
   >          .TRUE.. To select a complex conjugate pair of eigenvalues   
   >          w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,   
   >          either SELECT(j) or SELECT(j+1) or both must be set to   
   >          .TRUE.; a complex conjugate pair of eigenvalues must be   
   >          either both included in the cluster or both excluded.   
   > \endverbatim   
   >   
   > \param[in] N   
   > \verbatim   
   >          N is INTEGER   
   >          The order of the matrix T. N >= 0.   
   > \endverbatim   
   >   
   > \param[in,out] T   
   > \verbatim   
   >          T is DOUBLE PRECISION array, dimension (LDT,N)   
   >          On entry, the upper quasi-triangular matrix T, in Schur   
   >          canonical form.   
   >          On exit, T is overwritten by the reordered matrix T, again in   
   >          Schur canonical form, with the selected eigenvalues in the   
   >          leading diagonal blocks.   
   > \endverbatim   
   >   
   > \param[in] LDT   
   > \verbatim   
   >          LDT is INTEGER   
   >          The leading dimension of the array T. LDT >= max(1,N).   
   > \endverbatim   
   >   
   > \param[in,out] Q   
   > \verbatim   
   >          Q is DOUBLE PRECISION array, dimension (LDQ,N)   
   >          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.   
   >          On exit, if COMPQ = 'V', Q has been postmultiplied by the   
   >          orthogonal transformation matrix which reorders T; the   
   >          leading M columns of Q form an orthonormal basis for the   
   >          specified invariant subspace.   
   >          If COMPQ = 'N', Q is not referenced.   
   > \endverbatim   
   >   
   > \param[in] LDQ   
   > \verbatim   
   >          LDQ is INTEGER   
   >          The leading dimension of the array Q.   
   >          LDQ >= 1; and if COMPQ = 'V', LDQ >= N.   
   > \endverbatim   
   >   
   > \param[out] WR   
   > \verbatim   
   >          WR is DOUBLE PRECISION array, dimension (N)   
   > \endverbatim   
   > \param[out] WI   
   > \verbatim   
   >          WI is DOUBLE PRECISION array, dimension (N)   
   >   
   >          The real and imaginary parts, respectively, of the reordered   
   >          eigenvalues of T. The eigenvalues are stored in the same   
   >          order as on the diagonal of T, with WR(i) = T(i,i) and, if   
   >          T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and   
   >          WI(i+1) = -WI(i). Note that if a complex eigenvalue is   
   >          sufficiently ill-conditioned, then its value may differ   
   >          significantly from its value before reordering.   
   > \endverbatim   
   >   
   > \param[out] M   
   > \verbatim   
   >          M is INTEGER   
   >          The dimension of the specified invariant subspace.   
   >          0 < = M <= N.   
   > \endverbatim   
   >   
   > \param[out] S   
   > \verbatim   
   >          S is DOUBLE PRECISION   
   >          If JOB = 'E' or 'B', S is a lower bound on the reciprocal   
   >          condition number for the selected cluster of eigenvalues.   
   >          S cannot underestimate the true reciprocal condition number   
   >          by more than a factor of sqrt(N). If M = 0 or N, S = 1.   
   >          If JOB = 'N' or 'V', S is not referenced.   
   > \endverbatim   
   >   
   > \param[out] SEP   
   > \verbatim   
   >          SEP is DOUBLE PRECISION   
   >          If JOB = 'V' or 'B', SEP is the estimated reciprocal   
   >          condition number of the specified invariant subspace. If   
   >          M = 0 or N, SEP = norm(T).   
   >          If JOB = 'N' or 'E', SEP is not referenced.   
   > \endverbatim   
   >   
   > \param[out] WORK   
   > \verbatim   
   >          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))   
   >          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.   
   > \endverbatim   
   >   
   > \param[in] LWORK   
   > \verbatim   
   >          LWORK is INTEGER   
   >          The dimension of the array WORK.   
   >          If JOB = 'N', LWORK >= max(1,N);   
   >          if JOB = 'E', LWORK >= max(1,M*(N-M));   
   >          if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).   
   >   
   >          If LWORK = -1, then a workspace query is assumed; the routine   
   >          only calculates the optimal size of the WORK array, returns   
   >          this value as the first entry of the WORK array, and no error   
   >          message related to LWORK is issued by XERBLA.   
   > \endverbatim   
   >   
   > \param[out] IWORK   
   > \verbatim   
   >          IWORK is INTEGER array, dimension (MAX(1,LIWORK))   
   >          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.   
   > \endverbatim   
   >   
   > \param[in] LIWORK   
   > \verbatim   
   >          LIWORK is INTEGER   
   >          The dimension of the array IWORK.   
   >          If JOB = 'N' or 'E', LIWORK >= 1;   
   >          if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).   
   >   
   >          If LIWORK = -1, then a workspace query is assumed; the   
   >          routine only calculates the optimal size of the IWORK array,   
   >          returns this value as the first entry of the IWORK array, and   
   >          no error message related to LIWORK is issued by XERBLA.   
   > \endverbatim   
   >   
   > \param[out] INFO   
   > \verbatim   
   >          INFO is INTEGER   
   >          = 0: successful exit   
   >          < 0: if INFO = -i, the i-th argument had an illegal value   
   >          = 1: reordering of T failed because some eigenvalues are too   
   >               close to separate (the problem is very ill-conditioned);   
   >               T may have been partially reordered, and WR and WI   
   >               contain the eigenvalues in the same order as in T; S and   
   >               SEP (if requested) are set to zero.   
   > \endverbatim   

    Authors:   
    ========   

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

   > \date April 2012   

   > \ingroup doubleOTHERcomputational   

   > \par Further Details:   
    =====================   
   >   
   > \verbatim   
   >   
   >  DTRSEN first collects the selected eigenvalues by computing an   
   >  orthogonal transformation Z to move them to the top left corner of T.   
   >  In other words, the selected eigenvalues are the eigenvalues of T11   
   >  in:   
   >   
   >          Z**T * T * Z = ( T11 T12 ) n1   
   >                         (  0  T22 ) n2   
   >                            n1  n2   
   >   
   >  where N = n1+n2 and Z**T means the transpose of Z. The first n1 columns   
   >  of Z span the specified invariant subspace of T.   
   >   
   >  If T has been obtained from the real Schur factorization of a matrix   
   >  A = Q*T*Q**T, then the reordered real Schur factorization of A is given   
   >  by A = (Q*Z)*(Z**T*T*Z)*(Q*Z)**T, and the first n1 columns of Q*Z span   
   >  the corresponding invariant subspace of A.   
   >   
   >  The reciprocal condition number of the average of the eigenvalues of   
   >  T11 may be returned in S. S lies between 0 (very badly conditioned)   
   >  and 1 (very well conditioned). It is computed as follows. First we   
   >  compute R so that   
   >   
   >                         P = ( I  R ) n1   
   >                             ( 0  0 ) n2   
   >                               n1 n2   
   >   
   >  is the projector on the invariant subspace associated with T11.   
   >  R is the solution of the Sylvester equation:   
   >   
   >                        T11*R - R*T22 = T12.   
   >   
   >  Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote   
   >  the two-norm of M. Then S is computed as the lower bound   
   >   
   >                      (1 + F-norm(R)**2)**(-1/2)   
   >   
   >  on the reciprocal of 2-norm(P), the true reciprocal condition number.   
   >  S cannot underestimate 1 / 2-norm(P) by more than a factor of   
   >  sqrt(N).   
   >   
   >  An approximate error bound for the computed average of the   
   >  eigenvalues of T11 is   
   >   
   >                         EPS * norm(T) / S   
   >   
   >  where EPS is the machine precision.   
   >   
   >  The reciprocal condition number of the right invariant subspace   
   >  spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.   
   >  SEP is defined as the separation of T11 and T22:   
   >   
   >                     sep( T11, T22 ) = sigma-min( C )   
   >   
   >  where sigma-min(C) is the smallest singular value of the   
   >  n1*n2-by-n1*n2 matrix   
   >   
   >     C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )   
   >   
   >  I(m) is an m by m identity matrix, and kprod denotes the Kronecker   
   >  product. We estimate sigma-min(C) by the reciprocal of an estimate of   
   >  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)   
   >  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).   
   >   
   >  When SEP is small, small changes in T can cause large changes in   
   >  the invariant subspace. An approximate bound on the maximum angular   
   >  error in the computed right invariant subspace is   
   >   
   >                      EPS * norm(T) / SEP   
   > \endverbatim   
   >   
    =====================================================================   
   Subroutine */ int igraphdtrsen_(char *job, char *compq, logical *select, integer 
	*n, doublereal *t, integer *ldt, doublereal *q, integer *ldq, 
	doublereal *wr, doublereal *wi, integer *m, doublereal *s, doublereal 
	*sep, doublereal *work, integer *lwork, integer *iwork, integer *
	liwork, integer *info)
{
    /* System generated locals */
    integer q_dim1, q_offset, t_dim1, t_offset, i__1, i__2;
    doublereal d__1, d__2;

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

    /* Local variables */
    integer k, n1, n2, kk, nn, ks;
    doublereal est;
    integer kase;
    logical pair;
    integer ierr;
    logical swap;
    doublereal scale;
    extern logical igraphlsame_(char *, char *);
    integer isave[3], lwmin = 0;
    logical wantq, wants;
    doublereal rnorm;
    extern /* Subroutine */ int igraphdlacn2_(integer *, doublereal *, doublereal *,
	     integer *, doublereal *, integer *, integer *);
    extern doublereal igraphdlange_(char *, integer *, integer *, doublereal *, 
	    integer *, doublereal *);
    extern /* Subroutine */ int igraphdlacpy_(char *, integer *, integer *, 
	    doublereal *, integer *, doublereal *, integer *), 
	    igraphxerbla_(char *, integer *, ftnlen);
    logical wantbh;
    extern /* Subroutine */ int igraphdtrexc_(char *, integer *, doublereal *, 
	    integer *, doublereal *, integer *, integer *, integer *, 
	    doublereal *, integer *);
    integer liwmin;
    logical wantsp, lquery;
    extern /* Subroutine */ int igraphdtrsyl_(char *, char *, integer *, integer *, 
	    integer *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, integer *, doublereal *, integer *);


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


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


       Decode and test the input parameters   

       Parameter adjustments */
    --select;
    t_dim1 = *ldt;
    t_offset = 1 + t_dim1;
    t -= t_offset;
    q_dim1 = *ldq;
    q_offset = 1 + q_dim1;
    q -= q_offset;
    --wr;
    --wi;
    --work;
    --iwork;

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

    *info = 0;
    lquery = *lwork == -1;
    if (! igraphlsame_(job, "N") && ! wants && ! wantsp) {
	*info = -1;
    } else if (! igraphlsame_(compq, "N") && ! wantq) {
	*info = -2;
    } else if (*n < 0) {
	*info = -4;
    } else if (*ldt < max(1,*n)) {
	*info = -6;
    } else if (*ldq < 1 || wantq && *ldq < *n) {
	*info = -8;
    } else {

/*        Set M to the dimension of the specified invariant subspace,   
          and test LWORK and LIWORK. */

	*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: */
	}

	n1 = *m;
	n2 = *n - *m;
	nn = n1 * n2;

	if (wantsp) {
/* Computing MAX */
	    i__1 = 1, i__2 = nn << 1;
	    lwmin = max(i__1,i__2);
	    liwmin = max(1,nn);
	} else if (igraphlsame_(job, "N")) {
	    lwmin = max(1,*n);
	    liwmin = 1;
	} else if (igraphlsame_(job, "E")) {
	    lwmin = max(1,nn);
	    liwmin = 1;
	}

	if (*lwork < lwmin && ! lquery) {
	    *info = -15;
	} else if (*liwork < liwmin && ! lquery) {
	    *info = -17;
	}
    }

    if (*info == 0) {
	work[1] = (doublereal) lwmin;
	iwork[1] = liwmin;
    }

    if (*info != 0) {
	i__1 = -(*info);
	igraphxerbla_("DTRSEN", &i__1, (ftnlen)6);
	return 0;
    } else if (lquery) {
	return 0;
    }

/*     Quick return if possible. */

    if (*m == *n || *m == 0) {
	if (wants) {
	    *s = 1.;
	}
	if (wantsp) {
	    *sep = igraphdlange_("1", n, n, &t[t_offset], ldt, &work[1]);
	}
	goto L40;
    }

/*     Collect the selected blocks at the top-left corner of T. */

    ks = 0;
    pair = FALSE_;
    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {
	if (pair) {
	    pair = FALSE_;
	} else {
	    swap = select[k];
	    if (k < *n) {
		if (t[k + 1 + k * t_dim1] != 0.) {
		    pair = TRUE_;
		    swap = swap || select[k + 1];
		}
	    }
	    if (swap) {
		++ks;

/*              Swap the K-th block to position KS. */

		ierr = 0;
		kk = k;
		if (k != ks) {
		    igraphdtrexc_(compq, n, &t[t_offset], ldt, &q[q_offset], ldq, &
			    kk, &ks, &work[1], &ierr);
		}
		if (ierr == 1 || ierr == 2) {

/*                 Blocks too close to swap: exit. */

		    *info = 1;
		    if (wants) {
			*s = 0.;
		    }
		    if (wantsp) {
			*sep = 0.;
		    }
		    goto L40;
		}
		if (pair) {
		    ++ks;
		}
	    }
	}
/* L20: */
    }

    if (wants) {

/*        Solve Sylvester equation for R:   

             T11*R - R*T22 = scale*T12 */

	igraphdlacpy_("F", &n1, &n2, &t[(n1 + 1) * t_dim1 + 1], ldt, &work[1], &n1);
	igraphdtrsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 1 + (n1 
		+ 1) * t_dim1], ldt, &work[1], &n1, &scale, &ierr);

/*        Estimate the reciprocal of the condition number of the cluster   
          of eigenvalues. */

	rnorm = igraphdlange_("F", &n1, &n2, &work[1], &n1, &work[1]);
	if (rnorm == 0.) {
	    *s = 1.;
	} else {
	    *s = scale / (sqrt(scale * scale / rnorm + rnorm) * sqrt(rnorm));
	}
    }

    if (wantsp) {

/*        Estimate sep(T11,T22). */

	est = 0.;
	kase = 0;
L30:
	igraphdlacn2_(&nn, &work[nn + 1], &work[1], &iwork[1], &est, &kase, isave);
	if (kase != 0) {
	    if (kase == 1) {

/*              Solve  T11*R - R*T22 = scale*X. */

		igraphdtrsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 
			1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &
			ierr);
	    } else {

/*              Solve T11**T*R - R*T22**T = scale*X. */

		igraphdtrsyl_("T", "T", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 
			1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &
			ierr);
	    }
	    goto L30;
	}

	*sep = scale / est;
    }

L40:

/*     Store the output eigenvalues in WR and WI. */

    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {
	wr[k] = t[k + k * t_dim1];
	wi[k] = 0.;
/* L50: */
    }
    i__1 = *n - 1;
    for (k = 1; k <= i__1; ++k) {
	if (t[k + 1 + k * t_dim1] != 0.) {
	    wi[k] = sqrt((d__1 = t[k + (k + 1) * t_dim1], abs(d__1))) * sqrt((
		    d__2 = t[k + 1 + k * t_dim1], abs(d__2)));
	    wi[k + 1] = -wi[k];
	}
/* L60: */
    }

    work[1] = (doublereal) lwmin;
    iwork[1] = liwmin;

    return 0;

/*     End of DTRSEN */

} /* igraphdtrsen_ */