*> \brief \b ZLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method. * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download ZLAED0 + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, * IWORK, INFO ) * * .. Scalar Arguments .. * INTEGER INFO, LDQ, LDQS, N, QSIZ * .. * .. Array Arguments .. * INTEGER IWORK( * ) * DOUBLE PRECISION D( * ), E( * ), RWORK( * ) * COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> Using the divide and conquer method, ZLAED0 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. *> \endverbatim * * Arguments: * ========== * *> \param[in] QSIZ *> \verbatim *> QSIZ is INTEGER *> The dimension of the unitary matrix used to reduce *> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The dimension of the symmetric tridiagonal matrix. N >= 0. *> \endverbatim *> *> \param[in,out] D *> \verbatim *> D is DOUBLE PRECISION array, dimension (N) *> On entry, the diagonal elements of the tridiagonal matrix. *> On exit, the eigenvalues in ascending order. *> \endverbatim *> *> \param[in,out] E *> \verbatim *> E is DOUBLE PRECISION array, dimension (N-1) *> On entry, the off-diagonal elements of the tridiagonal matrix. *> On exit, E has been destroyed. *> \endverbatim *> *> \param[in,out] Q *> \verbatim *> Q is COMPLEX*16 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. *> \endverbatim *> *> \param[in] LDQ *> \verbatim *> LDQ is INTEGER *> The leading dimension of the array Q. LDQ >= max(1,N). *> \endverbatim *> *> \param[out] IWORK *> \verbatim *> IWORK is 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 ) *> \endverbatim *> *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, *> dimension (1 + 3*N + 2*N*lg N + 3*N**2) *> ( lg( N ) = smallest integer k *> such that 2^k >= N ) *> \endverbatim *> *> \param[out] QSTORE *> \verbatim *> QSTORE is COMPLEX*16 array, dimension (LDQS, N) *> Used to store parts of *> the eigenvector matrix when the updating matrix multiplies *> take place. *> \endverbatim *> *> \param[in] LDQS *> \verbatim *> LDQS is INTEGER *> The leading dimension of the array QSTORE. *> LDQS >= max(1,N). *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is 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). *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date December 2016 * *> \ingroup complex16OTHERcomputational * * ===================================================================== SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, $ IWORK, INFO ) * * -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * * .. Scalar Arguments .. INTEGER INFO, LDQ, LDQS, N, QSIZ * .. * .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION D( * ), E( * ), RWORK( * ) COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * ) * .. * * ===================================================================== * * Warning: N could be as big as QSIZ! * * .. Parameters .. DOUBLE PRECISION TWO PARAMETER ( TWO = 2.D+0 ) * .. * .. Local Scalars .. INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM, $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM, $ J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1, $ SPM1, SPM2, SUBMAT, SUBPBS, TLVLS DOUBLE PRECISION TEMP * .. * .. External Subroutines .. EXTERNAL DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7 * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, LOG, MAX * .. * .. Executable Statements .. * * Test the input parameters. * 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.LT.MAX( 0, N ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN INFO = -6 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZLAED0', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 ) * * Determine the size and placement of the submatrices, and save in * the leading elements of IWORK. * IWORK( 1 ) = N SUBPBS = 1 TLVLS = 0 10 CONTINUE IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN DO 20 J = SUBPBS, 1, -1 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 IWORK( 2*J-1 ) = IWORK( J ) / 2 20 CONTINUE TLVLS = TLVLS + 1 SUBPBS = 2*SUBPBS GO TO 10 END IF DO 30 J = 2, SUBPBS IWORK( J ) = IWORK( J ) + IWORK( J-1 ) 30 CONTINUE * * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1 * using rank-1 modifications (cuts). * SPM1 = SUBPBS - 1 DO 40 I = 1, SPM1 SUBMAT = IWORK( I ) + 1 SMM1 = SUBMAT - 1 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) ) D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) ) 40 CONTINUE * INDXQ = 4*N + 3 * * Set up workspaces for eigenvalues only/accumulate new vectors * routine * TEMP = LOG( DBLE( N ) ) / LOG( TWO ) LGN = INT( TEMP ) IF( 2**LGN.LT.N ) $ LGN = LGN + 1 IF( 2**LGN.LT.N ) $ LGN = LGN + 1 IPRMPT = INDXQ + N + 1 IPERM = IPRMPT + N*LGN IQPTR = IPERM + N*LGN IGIVPT = IQPTR + N + 2 IGIVCL = IGIVPT + N*LGN * IGIVNM = 1 IQ = IGIVNM + 2*N*LGN IWREM = IQ + N**2 + 1 * Initialize pointers DO 50 I = 0, SUBPBS IWORK( IPRMPT+I ) = 1 IWORK( IGIVPT+I ) = 1 50 CONTINUE IWORK( IQPTR ) = 1 * * Solve each submatrix eigenproblem at the bottom of the divide and * conquer tree. * CURR = 0 DO 70 I = 0, SPM1 IF( I.EQ.0 ) THEN SUBMAT = 1 MATSIZ = IWORK( 1 ) ELSE SUBMAT = IWORK( I ) + 1 MATSIZ = IWORK( I+1 ) - IWORK( I ) END IF LL = IQ - 1 + IWORK( IQPTR+CURR ) CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), $ RWORK( LL ), MATSIZ, RWORK, INFO ) CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ), $ MATSIZ, QSTORE( 1, SUBMAT ), LDQS, $ RWORK( IWREM ) ) IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2 CURR = CURR + 1 IF( INFO.GT.0 ) THEN INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 RETURN END IF K = 1 DO 60 J = SUBMAT, IWORK( I+1 ) IWORK( INDXQ+J ) = K K = K + 1 60 CONTINUE 70 CONTINUE * * Successively merge eigensystems of adjacent submatrices * into eigensystem for the corresponding larger matrix. * * while ( SUBPBS > 1 ) * CURLVL = 1 80 CONTINUE IF( SUBPBS.GT.1 ) THEN SPM2 = SUBPBS - 2 DO 90 I = 0, SPM2, 2 IF( I.EQ.0 ) THEN 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 = CURPRB + 1 END IF * * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2) * into an eigensystem of size MATSIZ. ZLAED7 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. * CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB, $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS, $ E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ), $ RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ), $ IWORK( IPERM ), IWORK( IGIVPT ), $ IWORK( IGIVCL ), RWORK( IGIVNM ), $ Q( 1, SUBMAT ), RWORK( IWREM ), $ IWORK( SUBPBS+1 ), INFO ) IF( INFO.GT.0 ) THEN INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 RETURN END IF IWORK( I / 2+1 ) = IWORK( I+2 ) 90 CONTINUE SUBPBS = SUBPBS / 2 CURLVL = CURLVL + 1 GO TO 80 END IF * * end while * * Re-merge the eigenvalues/vectors which were deflated at the final * merge step. * DO 100 I = 1, N J = IWORK( INDXQ+I ) RWORK( I ) = D( J ) CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 ) 100 CONTINUE CALL DCOPY( N, RWORK, 1, D, 1 ) * RETURN * * End of ZLAED0 * END