/* -- 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 doublereal c_b11 = 0.;
static doublereal c_b12 = 1.;
static integer c__12 = 12;
static integer c__2 = 2;
static integer c__49 = 49;
/* > \brief \b DHSEQR
=========== DOCUMENTATION ===========
Online html documentation available at
http://www.netlib.org/lapack/explore-html/
> \htmlonly
> Download DHSEQR + dependencies
> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhseqr.
f">
> [TGZ]</a>
> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhseqr.
f">
> [ZIP]</a>
> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhseqr.
f">
> [TXT]</a>
> \endhtmlonly
Definition:
===========
SUBROUTINE DHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
LDZ, WORK, LWORK, INFO )
INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N
CHARACTER COMPZ, JOB
DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ),
$ Z( LDZ, * )
> \par Purpose:
=============
>
> \verbatim
>
> DHSEQR computes the eigenvalues of a Hessenberg matrix H
> and, optionally, the matrices T and Z from the Schur decomposition
> H = Z T Z**T, where T is an upper quasi-triangular matrix (the
> Schur form), and Z is the orthogonal matrix of Schur vectors.
>
> Optionally Z may be postmultiplied into an input orthogonal
> matrix Q so that this routine can give the Schur factorization
> of a matrix A which has been reduced to the Hessenberg form H
> by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
> \endverbatim
Arguments:
==========
> \param[in] JOB
> \verbatim
> JOB is CHARACTER*1
> = 'E': compute eigenvalues only;
> = 'S': compute eigenvalues and the Schur form T.
> \endverbatim
>
> \param[in] COMPZ
> \verbatim
> COMPZ is CHARACTER*1
> = 'N': no Schur vectors are computed;
> = 'I': Z is initialized to the unit matrix and the matrix Z
> of Schur vectors of H is returned;
> = 'V': Z must contain an orthogonal matrix Q on entry, and
> the product Q*Z is returned.
> \endverbatim
>
> \param[in] N
> \verbatim
> N is INTEGER
> The order of the matrix H. N .GE. 0.
> \endverbatim
>
> \param[in] ILO
> \verbatim
> ILO is INTEGER
> \endverbatim
>
> \param[in] IHI
> \verbatim
> IHI is INTEGER
>
> It is assumed that H is already upper triangular in rows
> and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
> set by a previous call to DGEBAL, and then passed to ZGEHRD
> when the matrix output by DGEBAL is reduced to Hessenberg
> form. Otherwise ILO and IHI should be set to 1 and N
> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
> If N = 0, then ILO = 1 and IHI = 0.
> \endverbatim
>
> \param[in,out] H
> \verbatim
> H is DOUBLE PRECISION array, dimension (LDH,N)
> On entry, the upper Hessenberg matrix H.
> On exit, if INFO = 0 and JOB = 'S', then H contains the
> upper quasi-triangular matrix T from the Schur decomposition
> (the Schur form); 2-by-2 diagonal blocks (corresponding to
> complex conjugate pairs of eigenvalues) are returned in
> standard form, with H(i,i) = H(i+1,i+1) and
> H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
> contents of H are unspecified on exit. (The output value of
> H when INFO.GT.0 is given under the description of INFO
> below.)
>
> Unlike earlier versions of DHSEQR, this subroutine may
> explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
> or j = IHI+1, IHI+2, ... N.
> \endverbatim
>
> \param[in] LDH
> \verbatim
> LDH is INTEGER
> The leading dimension of the array H. LDH .GE. max(1,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 computed
> eigenvalues. If two eigenvalues are computed as a complex
> conjugate pair, they are stored in consecutive elements of
> WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
> WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
> the same order as on the diagonal of the Schur form returned
> in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
> diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
> WI(i+1) = -WI(i).
> \endverbatim
>
> \param[in,out] Z
> \verbatim
> Z is DOUBLE PRECISION array, dimension (LDZ,N)
> If COMPZ = 'N', Z is not referenced.
> If COMPZ = 'I', on entry Z need not be set and on exit,
> if INFO = 0, Z contains the orthogonal matrix Z of the Schur
> vectors of H. If COMPZ = 'V', on entry Z must contain an
> N-by-N matrix Q, which is assumed to be equal to the unit
> matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
> if INFO = 0, Z contains Q*Z.
> Normally Q is the orthogonal matrix generated by DORGHR
> after the call to DGEHRD which formed the Hessenberg matrix
> H. (The output value of Z when INFO.GT.0 is given under
> the description of INFO below.)
> \endverbatim
>
> \param[in] LDZ
> \verbatim
> LDZ is INTEGER
> The leading dimension of the array Z. if COMPZ = 'I' or
> COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1.
> \endverbatim
>
> \param[out] WORK
> \verbatim
> WORK is DOUBLE PRECISION array, dimension (LWORK)
> On exit, if INFO = 0, WORK(1) returns an estimate of
> the optimal value for LWORK.
> \endverbatim
>
> \param[in] LWORK
> \verbatim
> LWORK is INTEGER
> The dimension of the array WORK. LWORK .GE. max(1,N)
> is sufficient and delivers very good and sometimes
> optimal performance. However, LWORK as large as 11*N
> may be required for optimal performance. A workspace
> query is recommended to determine the optimal workspace
> size.
>
> If LWORK = -1, then DHSEQR does a workspace query.
> In this case, DHSEQR checks the input parameters and
> estimates the optimal workspace size for the given
> values of N, ILO and IHI. The estimate is returned
> in WORK(1). No error message related to LWORK is
> issued by XERBLA. Neither H nor Z are accessed.
> \endverbatim
>
> \param[out] INFO
> \verbatim
> INFO is INTEGER
> = 0: successful exit
> .LT. 0: if INFO = -i, the i-th argument had an illegal
> value
> .GT. 0: if INFO = i, DHSEQR failed to compute all of
> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
> and WI contain those eigenvalues which have been
> successfully computed. (Failures are rare.)
>
> If INFO .GT. 0 and JOB = 'E', then on exit, the
> remaining unconverged eigenvalues are the eigen-
> values of the upper Hessenberg matrix rows and
> columns ILO through INFO of the final, output
> value of H.
>
> If INFO .GT. 0 and JOB = 'S', then on exit
>
> (*) (initial value of H)*U = U*(final value of H)
>
> where U is an orthogonal matrix. The final
> value of H is upper Hessenberg and quasi-triangular
> in rows and columns INFO+1 through IHI.
>
> If INFO .GT. 0 and COMPZ = 'V', then on exit
>
> (final value of Z) = (initial value of Z)*U
>
> where U is the orthogonal matrix in (*) (regard-
> less of the value of JOB.)
>
> If INFO .GT. 0 and COMPZ = 'I', then on exit
> (final value of Z) = U
> where U is the orthogonal matrix in (*) (regard-
> less of the value of JOB.)
>
> If INFO .GT. 0 and COMPZ = 'N', then Z is not
> accessed.
> \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 Contributors:
==================
>
> Karen Braman and Ralph Byers, Department of Mathematics,
> University of Kansas, USA
> \par Further Details:
=====================
>
> \verbatim
>
> Default values supplied by
> ILAENV(ISPEC,'DHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
> It is suggested that these defaults be adjusted in order
> to attain best performance in each particular
> computational environment.
>
> ISPEC=12: The DLAHQR vs DLAQR0 crossover point.
> Default: 75. (Must be at least 11.)
>
> ISPEC=13: Recommended deflation window size.
> This depends on ILO, IHI and NS. NS is the
> number of simultaneous shifts returned
> by ILAENV(ISPEC=15). (See ISPEC=15 below.)
> The default for (IHI-ILO+1).LE.500 is NS.
> The default for (IHI-ILO+1).GT.500 is 3*NS/2.
>
> ISPEC=14: Nibble crossover point. (See IPARMQ for
> details.) Default: 14% of deflation window
> size.
>
> ISPEC=15: Number of simultaneous shifts in a multishift
> QR iteration.
>
> If IHI-ILO+1 is ...
>
> greater than ...but less ... the
> or equal to ... than default is
>
> 1 30 NS = 2(+)
> 30 60 NS = 4(+)
> 60 150 NS = 10(+)
> 150 590 NS = **
> 590 3000 NS = 64
> 3000 6000 NS = 128
> 6000 infinity NS = 256
>
> (+) By default some or all matrices of this order
> are passed to the implicit double shift routine
> DLAHQR and this parameter is ignored. See
> ISPEC=12 above and comments in IPARMQ for
> details.
>
> (**) The asterisks (**) indicate an ad-hoc
> function of N increasing from 10 to 64.
>
> ISPEC=16: Select structured matrix multiply.
> If the number of simultaneous shifts (specified
> by ISPEC=15) is less than 14, then the default
> for ISPEC=16 is 0. Otherwise the default for
> ISPEC=16 is 2.
> \endverbatim
> \par References:
================
>
> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
> 929--947, 2002.
> \n
> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
> of Matrix Analysis, volume 23, pages 948--973, 2002.
=====================================================================
Subroutine */ int igraphdhseqr_(char *job, char *compz, integer *n, integer *ilo,
integer *ihi, doublereal *h__, integer *ldh, doublereal *wr,
doublereal *wi, doublereal *z__, integer *ldz, doublereal *work,
integer *lwork, integer *info)
{
/* System generated locals */
address a__1[2];
integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2[2], i__3;
doublereal d__1;
char ch__1[2];
/* Builtin functions
Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
/* Local variables */
integer i__;
doublereal hl[2401] /* was [49][49] */;
integer kbot, nmin;
extern logical igraphlsame_(char *, char *);
logical initz;
doublereal workl[49];
logical wantt, wantz;
extern /* Subroutine */ int igraphdlaqr0_(logical *, logical *, integer *,
integer *, integer *, doublereal *, integer *, doublereal *,
doublereal *, integer *, integer *, doublereal *, integer *,
doublereal *, integer *, integer *), igraphdlahqr_(logical *, logical *,
integer *, integer *, integer *, doublereal *, integer *,
doublereal *, doublereal *, integer *, integer *, doublereal *,
integer *, integer *), igraphdlacpy_(char *, integer *, integer *,
doublereal *, integer *, doublereal *, integer *),
igraphdlaset_(char *, integer *, integer *, doublereal *, doublereal *,
doublereal *, integer *);
extern integer igraphilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
extern /* Subroutine */ int igraphxerbla_(char *, integer *, ftnlen);
logical lquery;
/* -- 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
=====================================================================
==== Matrices of order NTINY or smaller must be processed by
. DLAHQR because of insufficient subdiagonal scratch space.
. (This is a hard limit.) ====
==== NL allocates some local workspace to help small matrices
. through a rare DLAHQR failure. NL .GT. NTINY = 11 is
. required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
. mended. (The default value of NMIN is 75.) Using NL = 49
. allows up to six simultaneous shifts and a 16-by-16
. deflation window. ====
==== Decode and check the input parameters. ====
Parameter adjustments */
h_dim1 = *ldh;
h_offset = 1 + h_dim1;
h__ -= h_offset;
--wr;
--wi;
z_dim1 = *ldz;
z_offset = 1 + z_dim1;
z__ -= z_offset;
--work;
/* Function Body */
wantt = igraphlsame_(job, "S");
initz = igraphlsame_(compz, "I");
wantz = initz || igraphlsame_(compz, "V");
work[1] = (doublereal) max(1,*n);
lquery = *lwork == -1;
*info = 0;
if (! igraphlsame_(job, "E") && ! wantt) {
*info = -1;
} else if (! igraphlsame_(compz, "N") && ! wantz) {
*info = -2;
} else if (*n < 0) {
*info = -3;
} else if (*ilo < 1 || *ilo > max(1,*n)) {
*info = -4;
} else if (*ihi < min(*ilo,*n) || *ihi > *n) {
*info = -5;
} else if (*ldh < max(1,*n)) {
*info = -7;
} else if (*ldz < 1 || wantz && *ldz < max(1,*n)) {
*info = -11;
} else if (*lwork < max(1,*n) && ! lquery) {
*info = -13;
}
if (*info != 0) {
/* ==== Quick return in case of invalid argument. ==== */
i__1 = -(*info);
igraphxerbla_("DHSEQR", &i__1, (ftnlen)6);
return 0;
} else if (*n == 0) {
/* ==== Quick return in case N = 0; nothing to do. ==== */
return 0;
} else if (lquery) {
/* ==== Quick return in case of a workspace query ==== */
igraphdlaqr0_(&wantt, &wantz, n, ilo, ihi, &h__[h_offset], ldh, &wr[1], &wi[
1], ilo, ihi, &z__[z_offset], ldz, &work[1], lwork, info);
/* ==== Ensure reported workspace size is backward-compatible with
. previous LAPACK versions. ====
Computing MAX */
d__1 = (doublereal) max(1,*n);
work[1] = max(d__1,work[1]);
return 0;
} else {
/* ==== copy eigenvalues isolated by DGEBAL ==== */
i__1 = *ilo - 1;
for (i__ = 1; i__ <= i__1; ++i__) {
wr[i__] = h__[i__ + i__ * h_dim1];
wi[i__] = 0.;
/* L10: */
}
i__1 = *n;
for (i__ = *ihi + 1; i__ <= i__1; ++i__) {
wr[i__] = h__[i__ + i__ * h_dim1];
wi[i__] = 0.;
/* L20: */
}
/* ==== Initialize Z, if requested ==== */
if (initz) {
igraphdlaset_("A", n, n, &c_b11, &c_b12, &z__[z_offset], ldz)
;
}
/* ==== Quick return if possible ==== */
if (*ilo == *ihi) {
wr[*ilo] = h__[*ilo + *ilo * h_dim1];
wi[*ilo] = 0.;
return 0;
}
/* ==== DLAHQR/DLAQR0 crossover point ====
Writing concatenation */
i__2[0] = 1, a__1[0] = job;
i__2[1] = 1, a__1[1] = compz;
s_cat(ch__1, a__1, i__2, &c__2, (ftnlen)2);
nmin = igraphilaenv_(&c__12, "DHSEQR", ch__1, n, ilo, ihi, lwork, (ftnlen)6,
(ftnlen)2);
nmin = max(11,nmin);
/* ==== DLAQR0 for big matrices; DLAHQR for small ones ==== */
if (*n > nmin) {
igraphdlaqr0_(&wantt, &wantz, n, ilo, ihi, &h__[h_offset], ldh, &wr[1],
&wi[1], ilo, ihi, &z__[z_offset], ldz, &work[1], lwork,
info);
} else {
/* ==== Small matrix ==== */
igraphdlahqr_(&wantt, &wantz, n, ilo, ihi, &h__[h_offset], ldh, &wr[1],
&wi[1], ilo, ihi, &z__[z_offset], ldz, info);
if (*info > 0) {
/* ==== A rare DLAHQR failure! DLAQR0 sometimes succeeds
. when DLAHQR fails. ==== */
kbot = *info;
if (*n >= 49) {
/* ==== Larger matrices have enough subdiagonal scratch
. space to call DLAQR0 directly. ==== */
igraphdlaqr0_(&wantt, &wantz, n, ilo, &kbot, &h__[h_offset],
ldh, &wr[1], &wi[1], ilo, ihi, &z__[z_offset],
ldz, &work[1], lwork, info);
} else {
/* ==== Tiny matrices don't have enough subdiagonal
. scratch space to benefit from DLAQR0. Hence,
. tiny matrices must be copied into a larger
. array before calling DLAQR0. ==== */
igraphdlacpy_("A", n, n, &h__[h_offset], ldh, hl, &c__49);
hl[*n + 1 + *n * 49 - 50] = 0.;
i__1 = 49 - *n;
igraphdlaset_("A", &c__49, &i__1, &c_b11, &c_b11, &hl[(*n + 1) *
49 - 49], &c__49);
igraphdlaqr0_(&wantt, &wantz, &c__49, ilo, &kbot, hl, &c__49, &
wr[1], &wi[1], ilo, ihi, &z__[z_offset], ldz,
workl, &c__49, info);
if (wantt || *info != 0) {
igraphdlacpy_("A", n, n, hl, &c__49, &h__[h_offset], ldh);
}
}
}
}
/* ==== Clear out the trash, if necessary. ==== */
if ((wantt || *info != 0) && *n > 2) {
i__1 = *n - 2;
i__3 = *n - 2;
igraphdlaset_("L", &i__1, &i__3, &c_b11, &c_b11, &h__[h_dim1 + 3], ldh);
}
/* ==== Ensure reported workspace size is backward-compatible with
. previous LAPACK versions. ====
Computing MAX */
d__1 = (doublereal) max(1,*n);
work[1] = max(d__1,work[1]);
}
/* ==== End of DHSEQR ==== */
return 0;
} /* igraphdhseqr_ */