--- rpl/lapack/lapack/zlaed8.f 2010/01/26 15:22:45 1.1.1.1
+++ rpl/lapack/lapack/zlaed8.f 2014/01/27 09:28:37 1.14
@@ -1,11 +1,237 @@
+*> \brief \b ZLAED8 used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZLAED8 + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA,
+* Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR,
+* GIVCOL, GIVNUM, INFO )
+*
+* .. Scalar Arguments ..
+* INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
+* DOUBLE PRECISION RHO
+* ..
+* .. Array Arguments ..
+* INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
+* $ INDXQ( * ), PERM( * )
+* DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
+* $ Z( * )
+* COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZLAED8 merges the two sets of eigenvalues together into a single
+*> sorted set. Then it tries to deflate the size of the problem.
+*> There are two ways in which deflation can occur: when two or more
+*> eigenvalues are close together or if there is a tiny element in the
+*> Z vector. For each such occurrence the order of the related secular
+*> equation problem is reduced by one.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[out] K
+*> \verbatim
+*> K is INTEGER
+*> Contains the number of non-deflated eigenvalues.
+*> This is the order of the related secular equation.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The dimension of the symmetric tridiagonal matrix. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] QSIZ
+*> \verbatim
+*> QSIZ is INTEGER
+*> The dimension of the unitary matrix used to reduce
+*> the dense or band matrix to tridiagonal form.
+*> QSIZ >= N if ICOMPQ = 1.
+*> \endverbatim
+*>
+*> \param[in,out] Q
+*> \verbatim
+*> Q is COMPLEX*16 array, dimension (LDQ,N)
+*> On entry, Q contains the eigenvectors of the partially solved
+*> system which has been previously updated in matrix
+*> multiplies with other partially solved eigensystems.
+*> On exit, Q contains the trailing (N-K) updated eigenvectors
+*> (those which were deflated) in its last N-K columns.
+*> \endverbatim
+*>
+*> \param[in] LDQ
+*> \verbatim
+*> LDQ is INTEGER
+*> The leading dimension of the array Q. LDQ >= max( 1, N ).
+*> \endverbatim
+*>
+*> \param[in,out] D
+*> \verbatim
+*> D is DOUBLE PRECISION array, dimension (N)
+*> On entry, D contains the eigenvalues of the two submatrices to
+*> be combined. On exit, D contains the trailing (N-K) updated
+*> eigenvalues (those which were deflated) sorted into increasing
+*> order.
+*> \endverbatim
+*>
+*> \param[in,out] RHO
+*> \verbatim
+*> RHO is DOUBLE PRECISION
+*> Contains the off diagonal element associated with the rank-1
+*> cut which originally split the two submatrices which are now
+*> being recombined. RHO is modified during the computation to
+*> the value required by DLAED3.
+*> \endverbatim
+*>
+*> \param[in] CUTPNT
+*> \verbatim
+*> CUTPNT is INTEGER
+*> Contains the location of the last eigenvalue in the leading
+*> sub-matrix. MIN(1,N) <= CUTPNT <= N.
+*> \endverbatim
+*>
+*> \param[in] Z
+*> \verbatim
+*> Z is DOUBLE PRECISION array, dimension (N)
+*> On input this vector contains the updating vector (the last
+*> row of the first sub-eigenvector matrix and the first row of
+*> the second sub-eigenvector matrix). The contents of Z are
+*> destroyed during the updating process.
+*> \endverbatim
+*>
+*> \param[out] DLAMDA
+*> \verbatim
+*> DLAMDA is DOUBLE PRECISION array, dimension (N)
+*> Contains a copy of the first K eigenvalues which will be used
+*> by DLAED3 to form the secular equation.
+*> \endverbatim
+*>
+*> \param[out] Q2
+*> \verbatim
+*> Q2 is COMPLEX*16 array, dimension (LDQ2,N)
+*> If ICOMPQ = 0, Q2 is not referenced. Otherwise,
+*> Contains a copy of the first K eigenvectors which will be used
+*> by DLAED7 in a matrix multiply (DGEMM) to update the new
+*> eigenvectors.
+*> \endverbatim
+*>
+*> \param[in] LDQ2
+*> \verbatim
+*> LDQ2 is INTEGER
+*> The leading dimension of the array Q2. LDQ2 >= max( 1, N ).
+*> \endverbatim
+*>
+*> \param[out] W
+*> \verbatim
+*> W is DOUBLE PRECISION array, dimension (N)
+*> This will hold the first k values of the final
+*> deflation-altered z-vector and will be passed to DLAED3.
+*> \endverbatim
+*>
+*> \param[out] INDXP
+*> \verbatim
+*> INDXP is INTEGER array, dimension (N)
+*> This will contain the permutation used to place deflated
+*> values of D at the end of the array. On output INDXP(1:K)
+*> points to the nondeflated D-values and INDXP(K+1:N)
+*> points to the deflated eigenvalues.
+*> \endverbatim
+*>
+*> \param[out] INDX
+*> \verbatim
+*> INDX is INTEGER array, dimension (N)
+*> This will contain the permutation used to sort the contents of
+*> D into ascending order.
+*> \endverbatim
+*>
+*> \param[in] INDXQ
+*> \verbatim
+*> INDXQ is INTEGER array, dimension (N)
+*> This contains the permutation which separately sorts the two
+*> sub-problems in D into ascending order. Note that elements in
+*> the second half of this permutation must first have CUTPNT
+*> added to their values in order to be accurate.
+*> \endverbatim
+*>
+*> \param[out] PERM
+*> \verbatim
+*> PERM is INTEGER array, dimension (N)
+*> Contains the permutations (from deflation and sorting) to be
+*> applied to each eigenblock.
+*> \endverbatim
+*>
+*> \param[out] GIVPTR
+*> \verbatim
+*> GIVPTR is INTEGER
+*> Contains the number of Givens rotations which took place in
+*> this subproblem.
+*> \endverbatim
+*>
+*> \param[out] GIVCOL
+*> \verbatim
+*> GIVCOL is INTEGER array, dimension (2, N)
+*> Each pair of numbers indicates a pair of columns to take place
+*> in a Givens rotation.
+*> \endverbatim
+*>
+*> \param[out] GIVNUM
+*> \verbatim
+*> GIVNUM is DOUBLE PRECISION array, dimension (2, N)
+*> Each number indicates the S value to be used in the
+*> corresponding Givens rotation.
+*> \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 September 2012
+*
+*> \ingroup complex16OTHERcomputational
+*
+* =====================================================================
SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA,
$ Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR,
$ GIVCOL, GIVNUM, INFO )
*
-* -- LAPACK routine (version 3.2) --
+* -- LAPACK computational routine (version 3.4.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
+* September 2012
*
* .. Scalar Arguments ..
INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
@@ -19,116 +245,6 @@
COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * )
* ..
*
-* Purpose
-* =======
-*
-* ZLAED8 merges the two sets of eigenvalues together into a single
-* sorted set. Then it tries to deflate the size of the problem.
-* There are two ways in which deflation can occur: when two or more
-* eigenvalues are close together or if there is a tiny element in the
-* Z vector. For each such occurrence the order of the related secular
-* equation problem is reduced by one.
-*
-* Arguments
-* =========
-*
-* K (output) INTEGER
-* Contains the number of non-deflated eigenvalues.
-* This is the order of the related secular equation.
-*
-* N (input) INTEGER
-* The dimension of the symmetric tridiagonal matrix. N >= 0.
-*
-* QSIZ (input) INTEGER
-* The dimension of the unitary matrix used to reduce
-* the dense or band matrix to tridiagonal form.
-* QSIZ >= N if ICOMPQ = 1.
-*
-* Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
-* On entry, Q contains the eigenvectors of the partially solved
-* system which has been previously updated in matrix
-* multiplies with other partially solved eigensystems.
-* On exit, Q contains the trailing (N-K) updated eigenvectors
-* (those which were deflated) in its last N-K columns.
-*
-* LDQ (input) INTEGER
-* The leading dimension of the array Q. LDQ >= max( 1, N ).
-*
-* D (input/output) DOUBLE PRECISION array, dimension (N)
-* On entry, D contains the eigenvalues of the two submatrices to
-* be combined. On exit, D contains the trailing (N-K) updated
-* eigenvalues (those which were deflated) sorted into increasing
-* order.
-*
-* RHO (input/output) DOUBLE PRECISION
-* Contains the off diagonal element associated with the rank-1
-* cut which originally split the two submatrices which are now
-* being recombined. RHO is modified during the computation to
-* the value required by DLAED3.
-*
-* CUTPNT (input) INTEGER
-* Contains the location of the last eigenvalue in the leading
-* sub-matrix. MIN(1,N) <= CUTPNT <= N.
-*
-* Z (input) DOUBLE PRECISION array, dimension (N)
-* On input this vector contains the updating vector (the last
-* row of the first sub-eigenvector matrix and the first row of
-* the second sub-eigenvector matrix). The contents of Z are
-* destroyed during the updating process.
-*
-* DLAMDA (output) DOUBLE PRECISION array, dimension (N)
-* Contains a copy of the first K eigenvalues which will be used
-* by DLAED3 to form the secular equation.
-*
-* Q2 (output) COMPLEX*16 array, dimension (LDQ2,N)
-* If ICOMPQ = 0, Q2 is not referenced. Otherwise,
-* Contains a copy of the first K eigenvectors which will be used
-* by DLAED7 in a matrix multiply (DGEMM) to update the new
-* eigenvectors.
-*
-* LDQ2 (input) INTEGER
-* The leading dimension of the array Q2. LDQ2 >= max( 1, N ).
-*
-* W (output) DOUBLE PRECISION array, dimension (N)
-* This will hold the first k values of the final
-* deflation-altered z-vector and will be passed to DLAED3.
-*
-* INDXP (workspace) INTEGER array, dimension (N)
-* This will contain the permutation used to place deflated
-* values of D at the end of the array. On output INDXP(1:K)
-* points to the nondeflated D-values and INDXP(K+1:N)
-* points to the deflated eigenvalues.
-*
-* INDX (workspace) INTEGER array, dimension (N)
-* This will contain the permutation used to sort the contents of
-* D into ascending order.
-*
-* INDXQ (input) INTEGER array, dimension (N)
-* This contains the permutation which separately sorts the two
-* sub-problems in D into ascending order. Note that elements in
-* the second half of this permutation must first have CUTPNT
-* added to their values in order to be accurate.
-*
-* PERM (output) INTEGER array, dimension (N)
-* Contains the permutations (from deflation and sorting) to be
-* applied to each eigenblock.
-*
-* GIVPTR (output) INTEGER
-* Contains the number of Givens rotations which took place in
-* this subproblem.
-*
-* GIVCOL (output) INTEGER array, dimension (2, N)
-* Each pair of numbers indicates a pair of columns to take place
-* in a Givens rotation.
-*
-* GIVNUM (output) DOUBLE PRECISION array, dimension (2, N)
-* Each number indicates the S value to be used in the
-* corresponding Givens rotation.
-*
-* INFO (output) INTEGER
-* = 0: successful exit.
-* < 0: if INFO = -i, the i-th argument had an illegal value.
-*
* =====================================================================
*
* .. Parameters ..
@@ -174,6 +290,13 @@
RETURN
END IF
*
+* Need to initialize GIVPTR to O here in case of quick exit
+* to prevent an unspecified code behavior (usually sigfault)
+* when IWORK array on entry to *stedc is not zeroed
+* (or at least some IWORK entries which used in *laed7 for GIVPTR).
+*
+ GIVPTR = 0
+*
* Quick return if possible
*
IF( N.EQ.0 )
@@ -241,7 +364,6 @@
* components of Z are zero in this new basis.
*
K = 0
- GIVPTR = 0
K2 = N + 1
DO 60 J = 1, N
IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN