/src/igraph/vendor/lapack/dtrsen.c

Source (jump to first uncovered line)
/*  -- translated by f2c (version 20191129).
   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_ */


Coverage Report

Created: 2023-09-25 06:05

Line	Count	Source (jump to first uncovered line)
1		/* -- translated by f2c (version 20191129).
2		You must link the resulting object file with libf2c:
3		on Microsoft Windows system, link with libf2c.lib;
4		on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5		or, if you install libf2c.a in a standard place, with -lf2c -lm
6		-- in that order, at the end of the command line, as in
7		cc *.o -lf2c -lm
8		Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10		http://www.netlib.org/f2c/libf2c.zip
11		*/
12
13		#include "f2c.h"
14
15		/* Table of constant values */
16
17		static integer c_n1 = -1;
18
19		/* > \brief \b DTRSEN
20
21		=========== DOCUMENTATION ===========
22
23		Online html documentation available at
24		http://www.netlib.org/lapack/explore-html/
25
26		> \htmlonly
27		> Download DTRSEN + dependencies
28		> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsen.
29		f">
30		> [TGZ]</a>
31		> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsen.
32		f">
33		> [ZIP]</a>
34		> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsen.
35		f">
36		> [TXT]</a>
37		> \endhtmlonly
38
39		Definition:
40		===========
41
42		SUBROUTINE DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
43		M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
44
45		CHARACTER COMPQ, JOB
46		INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N
47		DOUBLE PRECISION S, SEP
48		LOGICAL SELECT( * )
49		INTEGER IWORK( * )
50		DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
51		$ WR( * )
52
53
54		> \par Purpose:
55		=============
56		>
57		> \verbatim
58		>
59		> DTRSEN reorders the real Schur factorization of a real matrix
60		> A = QTQ**T, so that a selected cluster of eigenvalues appears in
61		> the leading diagonal blocks of the upper quasi-triangular matrix T,
62		> and the leading columns of Q form an orthonormal basis of the
63		> corresponding right invariant subspace.
64		>
65		> Optionally the routine computes the reciprocal condition numbers of
66		> the cluster of eigenvalues and/or the invariant subspace.
67		>
68		> T must be in Schur canonical form (as returned by DHSEQR), that is,
69		> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
70		> 2-by-2 diagonal block has its diagonal elements equal and its
71		> off-diagonal elements of opposite sign.
72		> \endverbatim
73
74		Arguments:
75		==========
76
77		> \param[in] JOB
78		> \verbatim
79		> JOB is CHARACTER*1
80		> Specifies whether condition numbers are required for the
81		> cluster of eigenvalues (S) or the invariant subspace (SEP):
82		> = 'N': none;
83		> = 'E': for eigenvalues only (S);
84		> = 'V': for invariant subspace only (SEP);
85		> = 'B': for both eigenvalues and invariant subspace (S and
86		> SEP).
87		> \endverbatim
88		>
89		> \param[in] COMPQ
90		> \verbatim
91		> COMPQ is CHARACTER*1
92		> = 'V': update the matrix Q of Schur vectors;
93		> = 'N': do not update Q.
94		> \endverbatim
95		>
96		> \param[in] SELECT
97		> \verbatim
98		> SELECT is LOGICAL array, dimension (N)
99		> SELECT specifies the eigenvalues in the selected cluster. To
100		> select a real eigenvalue w(j), SELECT(j) must be set to
101		> .TRUE.. To select a complex conjugate pair of eigenvalues
102		> w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
103		> either SELECT(j) or SELECT(j+1) or both must be set to
104		> .TRUE.; a complex conjugate pair of eigenvalues must be
105		> either both included in the cluster or both excluded.
106		> \endverbatim
107		>
108		> \param[in] N
109		> \verbatim
110		> N is INTEGER
111		> The order of the matrix T. N >= 0.
112		> \endverbatim
113		>
114		> \param[in,out] T
115		> \verbatim
116		> T is DOUBLE PRECISION array, dimension (LDT,N)
117		> On entry, the upper quasi-triangular matrix T, in Schur
118		> canonical form.
119		> On exit, T is overwritten by the reordered matrix T, again in
120		> Schur canonical form, with the selected eigenvalues in the
121		> leading diagonal blocks.
122		> \endverbatim
123		>
124		> \param[in] LDT
125		> \verbatim
126		> LDT is INTEGER
127		> The leading dimension of the array T. LDT >= max(1,N).
128		> \endverbatim
129		>
130		> \param[in,out] Q
131		> \verbatim
132		> Q is DOUBLE PRECISION array, dimension (LDQ,N)
133		> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
134		> On exit, if COMPQ = 'V', Q has been postmultiplied by the
135		> orthogonal transformation matrix which reorders T; the
136		> leading M columns of Q form an orthonormal basis for the
137		> specified invariant subspace.
138		> If COMPQ = 'N', Q is not referenced.
139		> \endverbatim
140		>
141		> \param[in] LDQ
142		> \verbatim
143		> LDQ is INTEGER
144		> The leading dimension of the array Q.
145		> LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
146		> \endverbatim
147		>
148		> \param[out] WR
149		> \verbatim
150		> WR is DOUBLE PRECISION array, dimension (N)
151		> \endverbatim
152		> \param[out] WI
153		> \verbatim
154		> WI is DOUBLE PRECISION array, dimension (N)
155		>
156		> The real and imaginary parts, respectively, of the reordered
157		> eigenvalues of T. The eigenvalues are stored in the same
158		> order as on the diagonal of T, with WR(i) = T(i,i) and, if
159		> T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and
160		> WI(i+1) = -WI(i). Note that if a complex eigenvalue is
161		> sufficiently ill-conditioned, then its value may differ
162		> significantly from its value before reordering.
163		> \endverbatim
164		>
165		> \param[out] M
166		> \verbatim
167		> M is INTEGER
168		> The dimension of the specified invariant subspace.
169		> 0 < = M <= N.
170		> \endverbatim
171		>
172		> \param[out] S
173		> \verbatim
174		> S is DOUBLE PRECISION
175		> If JOB = 'E' or 'B', S is a lower bound on the reciprocal
176		> condition number for the selected cluster of eigenvalues.
177		> S cannot underestimate the true reciprocal condition number
178		> by more than a factor of sqrt(N). If M = 0 or N, S = 1.
179		> If JOB = 'N' or 'V', S is not referenced.
180		> \endverbatim
181		>
182		> \param[out] SEP
183		> \verbatim
184		> SEP is DOUBLE PRECISION
185		> If JOB = 'V' or 'B', SEP is the estimated reciprocal
186		> condition number of the specified invariant subspace. If
187		> M = 0 or N, SEP = norm(T).
188		> If JOB = 'N' or 'E', SEP is not referenced.
189		> \endverbatim
190		>
191		> \param[out] WORK
192		> \verbatim
193		> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
194		> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
195		> \endverbatim
196		>
197		> \param[in] LWORK
198		> \verbatim
199		> LWORK is INTEGER
200		> The dimension of the array WORK.
201		> If JOB = 'N', LWORK >= max(1,N);
202		> if JOB = 'E', LWORK >= max(1,M*(N-M));
203		> if JOB = 'V' or 'B', LWORK >= max(1,2M(N-M)).
204		>
205		> If LWORK = -1, then a workspace query is assumed; the routine
206		> only calculates the optimal size of the WORK array, returns
207		> this value as the first entry of the WORK array, and no error
208		> message related to LWORK is issued by XERBLA.
209		> \endverbatim
210		>
211		> \param[out] IWORK
212		> \verbatim
213		> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
214		> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
215		> \endverbatim
216		>
217		> \param[in] LIWORK
218		> \verbatim
219		> LIWORK is INTEGER
220		> The dimension of the array IWORK.
221		> If JOB = 'N' or 'E', LIWORK >= 1;
222		> if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).
223		>
224		> If LIWORK = -1, then a workspace query is assumed; the
225		> routine only calculates the optimal size of the IWORK array,
226		> returns this value as the first entry of the IWORK array, and
227		> no error message related to LIWORK is issued by XERBLA.
228		> \endverbatim
229		>
230		> \param[out] INFO
231		> \verbatim
232		> INFO is INTEGER
233		> = 0: successful exit
234		> < 0: if INFO = -i, the i-th argument had an illegal value
235		> = 1: reordering of T failed because some eigenvalues are too
236		> close to separate (the problem is very ill-conditioned);
237		> T may have been partially reordered, and WR and WI
238		> contain the eigenvalues in the same order as in T; S and
239		> SEP (if requested) are set to zero.
240		> \endverbatim
241
242		Authors:
243		========
244
245		> \author Univ. of Tennessee
246		> \author Univ. of California Berkeley
247		> \author Univ. of Colorado Denver
248		> \author NAG Ltd.
249
250		> \date April 2012
251
252		> \ingroup doubleOTHERcomputational
253
254		> \par Further Details:
255		=====================
256		>
257		> \verbatim
258		>
259		> DTRSEN first collects the selected eigenvalues by computing an
260		> orthogonal transformation Z to move them to the top left corner of T.
261		> In other words, the selected eigenvalues are the eigenvalues of T11
262		> in:
263		>
264		> Z*T T * Z = ( T11 T12 ) n1
265		> ( 0 T22 ) n2
266		> n1 n2
267		>
268		> where N = n1+n2 and Z**T means the transpose of Z. The first n1 columns
269		> of Z span the specified invariant subspace of T.
270		>
271		> If T has been obtained from the real Schur factorization of a matrix
272		> A = QTQ**T, then the reordered real Schur factorization of A is given
273		> by A = (QZ)(Z*TTZ)(QZ)T, and the first n1 columns of QZ span
274		> the corresponding invariant subspace of A.
275		>
276		> The reciprocal condition number of the average of the eigenvalues of
277		> T11 may be returned in S. S lies between 0 (very badly conditioned)
278		> and 1 (very well conditioned). It is computed as follows. First we
279		> compute R so that
280		>
281		> P = ( I R ) n1
282		> ( 0 0 ) n2
283		> n1 n2
284		>
285		> is the projector on the invariant subspace associated with T11.
286		> R is the solution of the Sylvester equation:
287		>
288		> T11R - RT22 = T12.
289		>
290		> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
291		> the two-norm of M. Then S is computed as the lower bound
292		>
293		> (1 + F-norm(R)2)(-1/2)
294		>
295		> on the reciprocal of 2-norm(P), the true reciprocal condition number.
296		> S cannot underestimate 1 / 2-norm(P) by more than a factor of
297		> sqrt(N).
298		>
299		> An approximate error bound for the computed average of the
300		> eigenvalues of T11 is
301		>
302		> EPS * norm(T) / S
303		>
304		> where EPS is the machine precision.
305		>
306		> The reciprocal condition number of the right invariant subspace
307		> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
308		> SEP is defined as the separation of T11 and T22:
309		>
310		> sep( T11, T22 ) = sigma-min( C )
311		>
312		> where sigma-min(C) is the smallest singular value of the
313		> n1n2-by-n1n2 matrix
314		>
315		> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
316		>
317		> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
318		> product. We estimate sigma-min(C) by the reciprocal of an estimate of
319		> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
320		> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
321		>
322		> When SEP is small, small changes in T can cause large changes in
323		> the invariant subspace. An approximate bound on the maximum angular
324		> error in the computed right invariant subspace is
325		>
326		> EPS * norm(T) / SEP
327		> \endverbatim
328		>
329		=====================================================================
330		Subroutine / int igraphdtrsen_(char job, char compq, logical select, integer
331		n, doublereal t, integer ldt, doublereal q, integer *ldq,
332		doublereal wr, doublereal wi, integer m, doublereal s, doublereal
333		sep, doublereal work, integer lwork, integer iwork, integer *
334		liwork, integer *info)
335	0	{
336		/* System generated locals */
337	0	integer q_dim1, q_offset, t_dim1, t_offset, i__1, i__2;
338	0	doublereal d__1, d__2;
339
340		/* Builtin functions */
341	0	double sqrt(doublereal);
342
343		/* Local variables */
344	0	integer k, n1, n2, kk, nn, ks;
345	0	doublereal est;
346	0	integer kase;
347	0	logical pair;
348	0	integer ierr;
349	0	logical swap;
350	0	doublereal scale;
351	0	extern logical igraphlsame_(char , char );
352	0	integer isave[3], lwmin = 0;
353	0	logical wantq, wants;
354	0	doublereal rnorm;
355	0	extern /* Subroutine / int igraphdlacn2_(integer , doublereal , doublereal ,
356	0	integer , doublereal , integer , integer );
357	0	extern doublereal igraphdlange_(char , integer , integer , doublereal ,
358	0	integer , doublereal );
359	0	extern /* Subroutine / int igraphdlacpy_(char , integer , integer ,
360	0	doublereal , integer , doublereal , integer ),
361	0	igraphxerbla_(char , integer , ftnlen);
362	0	logical wantbh;
363	0	extern /* Subroutine / int igraphdtrexc_(char , integer , doublereal ,
364	0	integer , doublereal , integer , integer , integer *,
365	0	doublereal , integer );
366	0	integer liwmin;
367	0	logical wantsp, lquery;
368	0	extern /* Subroutine / int igraphdtrsyl_(char , char , integer , integer *,
369	0	integer , doublereal , integer , doublereal , integer *,
370	0	doublereal , integer , doublereal , integer );
371
372
373		/* -- LAPACK computational routine (version 3.4.1) --
374		-- LAPACK is a software package provided by Univ. of Tennessee, --
375		-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
376		April 2012
377
378
379		=====================================================================
380
381
382		Decode and test the input parameters
383
384		Parameter adjustments */
385	0	--select;
386	0	t_dim1 = *ldt;
387	0	t_offset = 1 + t_dim1;
388	0	t -= t_offset;
389	0	q_dim1 = *ldq;
390	0	q_offset = 1 + q_dim1;
391	0	q -= q_offset;
392	0	--wr;
393	0	--wi;
394	0	--work;
395	0	--iwork;
396
397		/* Function Body */
398	0	wantbh = igraphlsame_(job, "B");
399	0	wants = igraphlsame_(job, "E") \|\| wantbh;
400	0	wantsp = igraphlsame_(job, "V") \|\| wantbh;
401	0	wantq = igraphlsame_(compq, "V");
402
403	0	*info = 0;
404	0	lquery = *lwork == -1;
405	0	if (! igraphlsame_(job, "N") && ! wants && ! wantsp) {
406	0	*info = -1;
407	0	} else if (! igraphlsame_(compq, "N") && ! wantq) {
408	0	*info = -2;
409	0	} else if (*n < 0) {
410	0	*info = -4;
411	0	} else if (ldt < max(1,n)) {
412	0	*info = -6;
413	0	} else if (ldq < 1 \|\| wantq && ldq < *n) {
414	0	*info = -8;
415	0	} else {
416
417		/* Set M to the dimension of the specified invariant subspace,
418		and test LWORK and LIWORK. */
419
420	0	*m = 0;
421	0	pair = FALSE_;
422	0	i__1 = *n;
423	0	for (k = 1; k <= i__1; ++k) {
424	0	if (pair) {
425	0	pair = FALSE_;
426	0	} else {
427	0	if (k < *n) {
428	0	if (t[k + 1 + k * t_dim1] == 0.) {
429	0	if (select[k]) {
430	0	++(*m);
431	0	}
432	0	} else {
433	0	pair = TRUE_;
434	0	if (select[k] \|\| select[k + 1]) {
435	0	*m += 2;
436	0	}
437	0	}
438	0	} else {
439	0	if (select[*n]) {
440	0	++(*m);
441	0	}
442	0	}
443	0	}
444		/* L10: */
445	0	}
446
447	0	n1 = *m;
448	0	n2 = n - m;
449	0	nn = n1 * n2;
450
451	0	if (wantsp) {
452		/* Computing MAX */
453	0	i__1 = 1, i__2 = nn << 1;
454	0	lwmin = max(i__1,i__2);
455	0	liwmin = max(1,nn);
456	0	} else if (igraphlsame_(job, "N")) {
457	0	lwmin = max(1,*n);
458	0	liwmin = 1;
459	0	} else if (igraphlsame_(job, "E")) {
460	0	lwmin = max(1,nn);
461	0	liwmin = 1;
462	0	}
463
464	0	if (*lwork < lwmin && ! lquery) {
465	0	*info = -15;
466	0	} else if (*liwork < liwmin && ! lquery) {
467	0	*info = -17;
468	0	}
469	0	}
470
471	0	if (*info == 0) {
472	0	work[1] = (doublereal) lwmin;
473	0	iwork[1] = liwmin;
474	0	}
475
476	0	if (*info != 0) {
477	0	i__1 = -(*info);
478	0	igraphxerbla_("DTRSEN", &i__1, (ftnlen)6);
479	0	return 0;
480	0	} else if (lquery) {
481	0	return 0;
482	0	}
483
484		/* Quick return if possible. */
485
486	0	if (m == n \|\| *m == 0) {
487	0	if (wants) {
488	0	*s = 1.;
489	0	}
490	0	if (wantsp) {
491	0	*sep = igraphdlange_("1", n, n, &t[t_offset], ldt, &work[1]);
492	0	}
493	0	goto L40;
494	0	}
495
496		/* Collect the selected blocks at the top-left corner of T. */
497
498	0	ks = 0;
499	0	pair = FALSE_;
500	0	i__1 = *n;
501	0	for (k = 1; k <= i__1; ++k) {
502	0	if (pair) {
503	0	pair = FALSE_;
504	0	} else {
505	0	swap = select[k];
506	0	if (k < *n) {
507	0	if (t[k + 1 + k * t_dim1] != 0.) {
508	0	pair = TRUE_;
509	0	swap = swap \|\| select[k + 1];
510	0	}
511	0	}
512	0	if (swap) {
513	0	++ks;
514
515		/* Swap the K-th block to position KS. */
516
517	0	ierr = 0;
518	0	kk = k;
519	0	if (k != ks) {
520	0	igraphdtrexc_(compq, n, &t[t_offset], ldt, &q[q_offset], ldq, &
521	0	kk, &ks, &work[1], &ierr);
522	0	}
523	0	if (ierr == 1 \|\| ierr == 2) {
524
525		/* Blocks too close to swap: exit. */
526
527	0	*info = 1;
528	0	if (wants) {
529	0	*s = 0.;
530	0	}
531	0	if (wantsp) {
532	0	*sep = 0.;
533	0	}
534	0	goto L40;
535	0	}
536	0	if (pair) {
537	0	++ks;
538	0	}
539	0	}
540	0	}
541		/* L20: */
542	0	}
543
544	0	if (wants) {
545
546		/* Solve Sylvester equation for R:
547
548		T11R - RT22 = scaleT12 /
549
550	0	igraphdlacpy_("F", &n1, &n2, &t[(n1 + 1) * t_dim1 + 1], ldt, &work[1], &n1);
551	0	igraphdtrsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 1 + (n1
552	0	+ 1) * t_dim1], ldt, &work[1], &n1, &scale, &ierr);
553
554		/* Estimate the reciprocal of the condition number of the cluster
555		of eigenvalues. */
556
557	0	rnorm = igraphdlange_("F", &n1, &n2, &work[1], &n1, &work[1]);
558	0	if (rnorm == 0.) {
559	0	*s = 1.;
560	0	} else {
561	0	s = scale / (sqrt(scale scale / rnorm + rnorm) * sqrt(rnorm));
562	0	}
563	0	}
564
565	0	if (wantsp) {
566
567		/* Estimate sep(T11,T22). */
568
569	0	est = 0.;
570	0	kase = 0;
571	0	L30:
572	0	igraphdlacn2_(&nn, &work[nn + 1], &work[1], &iwork[1], &est, &kase, isave);
573	0	if (kase != 0) {
574	0	if (kase == 1) {
575
576		/* Solve T11R - RT22 = scaleX. /
577
578	0	igraphdtrsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 +
579	0	1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &
580	0	ierr);
581	0	} else {
582
583		/* Solve T11*TR - RT22T = scaleX. */
584
585	0	igraphdtrsyl_("T", "T", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 +
586	0	1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &
587	0	ierr);
588	0	}
589	0	goto L30;
590	0	}
591
592	0	*sep = scale / est;
593	0	}
594
595	0	L40:
596
597		/* Store the output eigenvalues in WR and WI. */
598
599	0	i__1 = *n;
600	0	for (k = 1; k <= i__1; ++k) {
601	0	wr[k] = t[k + k * t_dim1];
602	0	wi[k] = 0.;
603		/* L50: */
604	0	}
605	0	i__1 = *n - 1;
606	0	for (k = 1; k <= i__1; ++k) {
607	0	if (t[k + 1 + k * t_dim1] != 0.) {
608	0	wi[k] = sqrt((d__1 = t[k + (k + 1) * t_dim1], abs(d__1))) * sqrt((
609	0	d__2 = t[k + 1 + k * t_dim1], abs(d__2)));
610	0	wi[k + 1] = -wi[k];
611	0	}
612		/* L60: */
613	0	}
614
615	0	work[1] = (doublereal) lwmin;
616	0	iwork[1] = liwmin;
617
618	0	return 0;
619
620		/* End of DTRSEN */
621
622	0	} /* igraphdtrsen_ */
623