--- rpl/lapack/lapack/zlaed8.f 2010/01/26 15:22:45 1.1.1.1 +++ rpl/lapack/lapack/zlaed8.f 2023/08/07 08:39:28 1.19 @@ -1,11 +1,234 @@ +*> \brief \b ZLAED8 used by ZSTEDC. 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. +* +*> \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 -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 * * .. Scalar Arguments .. INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ @@ -19,116 +242,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 +287,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 +361,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