#include "blaswrap.h"
/* claed0.f -- translated by f2c (version 20061008).
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__9 = 9;
static integer c__0 = 0;
static integer c__2 = 2;
static integer c__1 = 1;
/* Subroutine */ int claed0_(integer *qsiz, integer *n, real *d__, real *e,
complex *q, integer *ldq, complex *qstore, integer *ldqs, real *rwork,
integer *iwork, integer *info)
{
/* System generated locals */
integer q_dim1, q_offset, qstore_dim1, qstore_offset, i__1, i__2;
real r__1;
/* Builtin functions */
double log(doublereal);
integer pow_ii(integer *, integer *);
/* Local variables */
static integer i__, j, k, ll, iq, lgn, msd2, smm1, spm1, spm2;
static real temp;
static integer curr, iperm;
extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
complex *, integer *);
static integer indxq, iwrem;
extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
integer *);
static integer iqptr;
extern /* Subroutine */ int claed7_(integer *, integer *, integer *,
integer *, integer *, integer *, real *, complex *, integer *,
real *, integer *, real *, integer *, integer *, integer *,
integer *, integer *, real *, complex *, real *, integer *,
integer *);
static integer tlvls;
extern /* Subroutine */ int clacrm_(integer *, integer *, complex *,
integer *, real *, integer *, complex *, integer *, real *);
static integer igivcl;
extern /* Subroutine */ int xerbla_(char *, integer *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
static integer igivnm, submat, curprb, subpbs, igivpt, curlvl, matsiz,
iprmpt, smlsiz;
extern /* Subroutine */ int ssteqr_(char *, integer *, real *, real *,
real *, integer *, real *, integer *);
/* -- LAPACK routine (version 3.1) --
Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
November 2006
Purpose
=======
Using the divide and conquer method, CLAED0 computes all eigenvalues
of a symmetric tridiagonal matrix which is one diagonal block of
those from reducing a dense or band Hermitian matrix and
corresponding eigenvectors of the dense or band matrix.
Arguments
=========
QSIZ (input) INTEGER
The dimension of the unitary matrix used to reduce
the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
N (input) INTEGER
The dimension of the symmetric tridiagonal matrix. N >= 0.
D (input/output) REAL array, dimension (N)
On entry, the diagonal elements of the tridiagonal matrix.
On exit, the eigenvalues in ascending order.
E (input/output) REAL array, dimension (N-1)
On entry, the off-diagonal elements of the tridiagonal matrix.
On exit, E has been destroyed.
Q (input/output) COMPLEX array, dimension (LDQ,N)
On entry, Q must contain an QSIZ x N matrix whose columns
unitarily orthonormal. It is a part of the unitary matrix
that reduces the full dense Hermitian matrix to a
(reducible) symmetric tridiagonal matrix.
LDQ (input) INTEGER
The leading dimension of the array Q. LDQ >= max(1,N).
IWORK (workspace) INTEGER array,
the dimension of IWORK must be at least
6 + 6*N + 5*N*lg N
( lg( N ) = smallest integer k
such that 2^k >= N )
RWORK (workspace) REAL array,
dimension (1 + 3*N + 2*N*lg N + 3*N**2)
( lg( N ) = smallest integer k
such that 2^k >= N )
QSTORE (workspace) COMPLEX array, dimension (LDQS, N)
Used to store parts of
the eigenvector matrix when the updating matrix multiplies
take place.
LDQS (input) INTEGER
The leading dimension of the array QSTORE.
LDQS >= max(1,N).
INFO (output) INTEGER
= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an illegal value.
> 0: The algorithm failed to compute an eigenvalue while
working on the submatrix lying in rows and columns
INFO/(N+1) through mod(INFO,N+1).
=====================================================================
Warning: N could be as big as QSIZ!
Test the input parameters.
Parameter adjustments */
--d__;
--e;
q_dim1 = *ldq;
q_offset = 1 + q_dim1;
q -= q_offset;
qstore_dim1 = *ldqs;
qstore_offset = 1 + qstore_dim1;
qstore -= qstore_offset;
--rwork;
--iwork;
/* Function Body */
*info = 0;
/* IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
INFO = -1
ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
$ THEN */
if (*qsiz < max(0,*n)) {
*info = -1;
} else if (*n < 0) {
*info = -2;
} else if (*ldq < max(1,*n)) {
*info = -6;
} else if (*ldqs < max(1,*n)) {
*info = -8;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("CLAED0", &i__1);
return 0;
}
/* Quick return if possible */
if (*n == 0) {
return 0;
}
smlsiz = ilaenv_(&c__9, "CLAED0", " ", &c__0, &c__0, &c__0, &c__0, (
ftnlen)6, (ftnlen)1);
/* Determine the size and placement of the submatrices, and save in
the leading elements of IWORK. */
iwork[1] = *n;
subpbs = 1;
tlvls = 0;
L10:
if (iwork[subpbs] > smlsiz) {
for (j = subpbs; j >= 1; --j) {
iwork[j * 2] = (iwork[j] + 1) / 2;
iwork[(j << 1) - 1] = iwork[j] / 2;
/* L20: */
}
++tlvls;
subpbs <<= 1;
goto L10;
}
i__1 = subpbs;
for (j = 2; j <= i__1; ++j) {
iwork[j] += iwork[j - 1];
/* L30: */
}
/* Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
using rank-1 modifications (cuts). */
spm1 = subpbs - 1;
i__1 = spm1;
for (i__ = 1; i__ <= i__1; ++i__) {
submat = iwork[i__] + 1;
smm1 = submat - 1;
d__[smm1] -= (r__1 = e[smm1], dabs(r__1));
d__[submat] -= (r__1 = e[smm1], dabs(r__1));
/* L40: */
}
indxq = (*n << 2) + 3;
/* Set up workspaces for eigenvalues only/accumulate new vectors
routine */
temp = log((real) (*n)) / log(2.f);
lgn = (integer) temp;
if (pow_ii(&c__2, &lgn) < *n) {
++lgn;
}
if (pow_ii(&c__2, &lgn) < *n) {
++lgn;
}
iprmpt = indxq + *n + 1;
iperm = iprmpt + *n * lgn;
iqptr = iperm + *n * lgn;
igivpt = iqptr + *n + 2;
igivcl = igivpt + *n * lgn;
igivnm = 1;
iq = igivnm + (*n << 1) * lgn;
/* Computing 2nd power */
i__1 = *n;
iwrem = iq + i__1 * i__1 + 1;
/* Initialize pointers */
i__1 = subpbs;
for (i__ = 0; i__ <= i__1; ++i__) {
iwork[iprmpt + i__] = 1;
iwork[igivpt + i__] = 1;
/* L50: */
}
iwork[iqptr] = 1;
/* Solve each submatrix eigenproblem at the bottom of the divide and
conquer tree. */
curr = 0;
i__1 = spm1;
for (i__ = 0; i__ <= i__1; ++i__) {
if (i__ == 0) {
submat = 1;
matsiz = iwork[1];
} else {
submat = iwork[i__] + 1;
matsiz = iwork[i__ + 1] - iwork[i__];
}
ll = iq - 1 + iwork[iqptr + curr];
ssteqr_("I", &matsiz, &d__[submat], &e[submat], &rwork[ll], &matsiz, &
rwork[1], info);
clacrm_(qsiz, &matsiz, &q[submat * q_dim1 + 1], ldq, &rwork[ll], &
matsiz, &qstore[submat * qstore_dim1 + 1], ldqs, &rwork[iwrem] );
/* Computing 2nd power */
i__2 = matsiz;
iwork[iqptr + curr + 1] = iwork[iqptr + curr] + i__2 * i__2;
++curr;
if (*info > 0) {
*info = submat * (*n + 1) + submat + matsiz - 1;
return 0;
}
k = 1;
i__2 = iwork[i__ + 1];
for (j = submat; j <= i__2; ++j) {
iwork[indxq + j] = k;
++k;
/* L60: */
}
/* L70: */
}
/* Successively merge eigensystems of adjacent submatrices
into eigensystem for the corresponding larger matrix.
while ( SUBPBS > 1 ) */
curlvl = 1;
L80:
if (subpbs > 1) {
spm2 = subpbs - 2;
i__1 = spm2;
for (i__ = 0; i__ <= i__1; i__ += 2) {
if (i__ == 0) {
submat = 1;
matsiz = iwork[2];
msd2 = iwork[1];
curprb = 0;
} else {
submat = iwork[i__] + 1;
matsiz = iwork[i__ + 2] - iwork[i__];
msd2 = matsiz / 2;
++curprb;
}
/* Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
into an eigensystem of size MATSIZ. CLAED7 handles the case
when the eigenvectors of a full or band Hermitian matrix (which
was reduced to tridiagonal form) are desired.
I am free to use Q as a valuable working space until Loop 150. */
claed7_(&matsiz, &msd2, qsiz, &tlvls, &curlvl, &curprb, &d__[
submat], &qstore[submat * qstore_dim1 + 1], ldqs, &e[
submat + msd2 - 1], &iwork[indxq + submat], &rwork[iq], &
iwork[iqptr], &iwork[iprmpt], &iwork[iperm], &iwork[
igivpt], &iwork[igivcl], &rwork[igivnm], &q[submat *
q_dim1 + 1], &rwork[iwrem], &iwork[subpbs + 1], info);
if (*info > 0) {
*info = submat * (*n + 1) + submat + matsiz - 1;
return 0;
}
iwork[i__ / 2 + 1] = iwork[i__ + 2];
/* L90: */
}
subpbs /= 2;
++curlvl;
goto L80;
}
/* end while
Re-merge the eigenvalues/vectors which were deflated at the final
merge step. */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
j = iwork[indxq + i__];
rwork[i__] = d__[j];
ccopy_(qsiz, &qstore[j * qstore_dim1 + 1], &c__1, &q[i__ * q_dim1 + 1]
, &c__1);
/* L100: */
}
scopy_(n, &rwork[1], &c__1, &d__[1], &c__1);
return 0;
/* End of CLAED0 */
} /* claed0_ */