--- rpl/lapack/lapack/zstemr.f 2010/12/21 13:53:55 1.7 +++ rpl/lapack/lapack/zstemr.f 2011/11/21 20:43:21 1.8 @@ -1,14 +1,338 @@ +*> \brief \b ZSTEMR +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZSTEMR + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, +* M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, +* IWORK, LIWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBZ, RANGE +* LOGICAL TRYRAC +* INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N +* DOUBLE PRECISION VL, VU +* .. +* .. Array Arguments .. +* INTEGER ISUPPZ( * ), IWORK( * ) +* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) +* COMPLEX*16 Z( LDZ, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZSTEMR computes selected eigenvalues and, optionally, eigenvectors +*> of a real symmetric tridiagonal matrix T. Any such unreduced matrix has +*> a well defined set of pairwise different real eigenvalues, the corresponding +*> real eigenvectors are pairwise orthogonal. +*> +*> The spectrum may be computed either completely or partially by specifying +*> either an interval (VL,VU] or a range of indices IL:IU for the desired +*> eigenvalues. +*> +*> Depending on the number of desired eigenvalues, these are computed either +*> by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are +*> computed by the use of various suitable L D L^T factorizations near clusters +*> of close eigenvalues (referred to as RRRs, Relatively Robust +*> Representations). An informal sketch of the algorithm follows. +*> +*> For each unreduced block (submatrix) of T, +*> (a) Compute T - sigma I = L D L^T, so that L and D +*> define all the wanted eigenvalues to high relative accuracy. +*> This means that small relative changes in the entries of D and L +*> cause only small relative changes in the eigenvalues and +*> eigenvectors. The standard (unfactored) representation of the +*> tridiagonal matrix T does not have this property in general. +*> (b) Compute the eigenvalues to suitable accuracy. +*> If the eigenvectors are desired, the algorithm attains full +*> accuracy of the computed eigenvalues only right before +*> the corresponding vectors have to be computed, see steps c) and d). +*> (c) For each cluster of close eigenvalues, select a new +*> shift close to the cluster, find a new factorization, and refine +*> the shifted eigenvalues to suitable accuracy. +*> (d) For each eigenvalue with a large enough relative separation compute +*> the corresponding eigenvector by forming a rank revealing twisted +*> factorization. Go back to (c) for any clusters that remain. +*> +*> For more details, see: +*> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations +*> to compute orthogonal eigenvectors of symmetric tridiagonal matrices," +*> Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. +*> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and +*> Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, +*> 2004. Also LAPACK Working Note 154. +*> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric +*> tridiagonal eigenvalue/eigenvector problem", +*> Computer Science Division Technical Report No. UCB/CSD-97-971, +*> UC Berkeley, May 1997. +*> +*> Further Details +*> 1.ZSTEMR works only on machines which follow IEEE-754 +*> floating-point standard in their handling of infinities and NaNs. +*> This permits the use of efficient inner loops avoiding a check for +*> zero divisors. +*> +*> 2. LAPACK routines can be used to reduce a complex Hermitean matrix to +*> real symmetric tridiagonal form. +*> +*> (Any complex Hermitean tridiagonal matrix has real values on its diagonal +*> and potentially complex numbers on its off-diagonals. By applying a +*> similarity transform with an appropriate diagonal matrix +*> diag(1,e^{i \phy_1}, ... , e^{i \phy_{n-1}}), the complex Hermitean +*> matrix can be transformed into a real symmetric matrix and complex +*> arithmetic can be entirely avoided.) +*> +*> While the eigenvectors of the real symmetric tridiagonal matrix are real, +*> the eigenvectors of original complex Hermitean matrix have complex entries +*> in general. +*> Since LAPACK drivers overwrite the matrix data with the eigenvectors, +*> ZSTEMR accepts complex workspace to facilitate interoperability +*> with ZUNMTR or ZUPMTR. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBZ +*> \verbatim +*> JOBZ is CHARACTER*1 +*> = 'N': Compute eigenvalues only; +*> = 'V': Compute eigenvalues and eigenvectors. +*> \endverbatim +*> +*> \param[in] RANGE +*> \verbatim +*> RANGE is CHARACTER*1 +*> = 'A': all eigenvalues will be found. +*> = 'V': all eigenvalues in the half-open interval (VL,VU] +*> will be found. +*> = 'I': the IL-th through IU-th eigenvalues will be found. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix. N >= 0. +*> \endverbatim +*> +*> \param[in,out] D +*> \verbatim +*> D is DOUBLE PRECISION array, dimension (N) +*> On entry, the N diagonal elements of the tridiagonal matrix +*> T. On exit, D is overwritten. +*> \endverbatim +*> +*> \param[in,out] E +*> \verbatim +*> E is DOUBLE PRECISION array, dimension (N) +*> On entry, the (N-1) subdiagonal elements of the tridiagonal +*> matrix T in elements 1 to N-1 of E. E(N) need not be set on +*> input, but is used internally as workspace. +*> On exit, E is overwritten. +*> \endverbatim +*> +*> \param[in] VL +*> \verbatim +*> VL is DOUBLE PRECISION +*> \endverbatim +*> +*> \param[in] VU +*> \verbatim +*> VU is DOUBLE PRECISION +*> +*> If RANGE='V', the lower and upper bounds of the interval to +*> be searched for eigenvalues. VL < VU. +*> Not referenced if RANGE = 'A' or 'I'. +*> \endverbatim +*> +*> \param[in] IL +*> \verbatim +*> IL is INTEGER +*> \endverbatim +*> +*> \param[in] IU +*> \verbatim +*> IU is INTEGER +*> +*> If RANGE='I', the indices (in ascending order) of the +*> smallest and largest eigenvalues to be returned. +*> 1 <= IL <= IU <= N, if N > 0. +*> Not referenced if RANGE = 'A' or 'V'. +*> \endverbatim +*> +*> \param[out] M +*> \verbatim +*> M is INTEGER +*> The total number of eigenvalues found. 0 <= M <= N. +*> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. +*> \endverbatim +*> +*> \param[out] W +*> \verbatim +*> W is DOUBLE PRECISION array, dimension (N) +*> The first M elements contain the selected eigenvalues in +*> ascending order. +*> \endverbatim +*> +*> \param[out] Z +*> \verbatim +*> Z is COMPLEX*16 array, dimension (LDZ, max(1,M) ) +*> If JOBZ = 'V', and if INFO = 0, then the first M columns of Z +*> contain the orthonormal eigenvectors of the matrix T +*> corresponding to the selected eigenvalues, with the i-th +*> column of Z holding the eigenvector associated with W(i). +*> If JOBZ = 'N', then Z is not referenced. +*> Note: the user must ensure that at least max(1,M) columns are +*> supplied in the array Z; if RANGE = 'V', the exact value of M +*> is not known in advance and can be computed with a workspace +*> query by setting NZC = -1, see below. +*> \endverbatim +*> +*> \param[in] LDZ +*> \verbatim +*> LDZ is INTEGER +*> The leading dimension of the array Z. LDZ >= 1, and if +*> JOBZ = 'V', then LDZ >= max(1,N). +*> \endverbatim +*> +*> \param[in] NZC +*> \verbatim +*> NZC is INTEGER +*> The number of eigenvectors to be held in the array Z. +*> If RANGE = 'A', then NZC >= max(1,N). +*> If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU]. +*> If RANGE = 'I', then NZC >= IU-IL+1. +*> If NZC = -1, then a workspace query is assumed; the +*> routine calculates the number of columns of the array Z that +*> are needed to hold the eigenvectors. +*> This value is returned as the first entry of the Z array, and +*> no error message related to NZC is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] ISUPPZ +*> \verbatim +*> ISUPPZ is INTEGER ARRAY, dimension ( 2*max(1,M) ) +*> The support of the eigenvectors in Z, i.e., the indices +*> indicating the nonzero elements in Z. The i-th computed eigenvector +*> is nonzero only in elements ISUPPZ( 2*i-1 ) through +*> ISUPPZ( 2*i ). This is relevant in the case when the matrix +*> is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0. +*> \endverbatim +*> +*> \param[in,out] TRYRAC +*> \verbatim +*> TRYRAC is LOGICAL +*> If TRYRAC.EQ..TRUE., indicates that the code should check whether +*> the tridiagonal matrix defines its eigenvalues to high relative +*> accuracy. If so, the code uses relative-accuracy preserving +*> algorithms that might be (a bit) slower depending on the matrix. +*> If the matrix does not define its eigenvalues to high relative +*> accuracy, the code can uses possibly faster algorithms. +*> If TRYRAC.EQ..FALSE., the code is not required to guarantee +*> relatively accurate eigenvalues and can use the fastest possible +*> techniques. +*> On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix +*> does not define its eigenvalues to high relative accuracy. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (LWORK) +*> On exit, if INFO = 0, WORK(1) returns the optimal +*> (and minimal) LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK. LWORK >= max(1,18*N) +*> if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK array, returns +*> this value as the first entry of the WORK array, and no error +*> message related to LWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (LIWORK) +*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. +*> \endverbatim +*> +*> \param[in] LIWORK +*> \verbatim +*> LIWORK is INTEGER +*> The dimension of the array IWORK. LIWORK >= max(1,10*N) +*> if the eigenvectors are desired, and LIWORK >= max(1,8*N) +*> if only the eigenvalues are to be computed. +*> If LIWORK = -1, then a workspace query is assumed; the +*> routine only calculates the optimal size of the IWORK array, +*> returns this value as the first entry of the IWORK array, and +*> no error message related to LIWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> On exit, INFO +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = 1X, internal error in DLARRE, +*> if INFO = 2X, internal error in ZLARRV. +*> Here, the digit X = ABS( IINFO ) < 10, where IINFO is +*> the nonzero error code returned by DLARRE or +*> ZLARRV, respectively. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complex16OTHERcomputational +* +*> \par Contributors: +* ================== +*> +*> Beresford Parlett, University of California, Berkeley, USA \n +*> Jim Demmel, University of California, Berkeley, USA \n +*> Inderjit Dhillon, University of Texas, Austin, USA \n +*> Osni Marques, LBNL/NERSC, USA \n +*> Christof Voemel, University of California, Berkeley, USA \n +* +* ===================================================================== SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, $ M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, $ IWORK, LIWORK, INFO ) - IMPLICIT NONE -* -* -- LAPACK computational routine (version 3.2.1) -- -* -* -- April 2009 -- * +* -- 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 * * .. Scalar Arguments .. CHARACTER JOBZ, RANGE @@ -22,215 +346,6 @@ COMPLEX*16 Z( LDZ, * ) * .. * -* Purpose -* ======= -* -* ZSTEMR computes selected eigenvalues and, optionally, eigenvectors -* of a real symmetric tridiagonal matrix T. Any such unreduced matrix has -* a well defined set of pairwise different real eigenvalues, the corresponding -* real eigenvectors are pairwise orthogonal. -* -* The spectrum may be computed either completely or partially by specifying -* either an interval (VL,VU] or a range of indices IL:IU for the desired -* eigenvalues. -* -* Depending on the number of desired eigenvalues, these are computed either -* by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are -* computed by the use of various suitable L D L^T factorizations near clusters -* of close eigenvalues (referred to as RRRs, Relatively Robust -* Representations). An informal sketch of the algorithm follows. -* -* For each unreduced block (submatrix) of T, -* (a) Compute T - sigma I = L D L^T, so that L and D -* define all the wanted eigenvalues to high relative accuracy. -* This means that small relative changes in the entries of D and L -* cause only small relative changes in the eigenvalues and -* eigenvectors. The standard (unfactored) representation of the -* tridiagonal matrix T does not have this property in general. -* (b) Compute the eigenvalues to suitable accuracy. -* If the eigenvectors are desired, the algorithm attains full -* accuracy of the computed eigenvalues only right before -* the corresponding vectors have to be computed, see steps c) and d). -* (c) For each cluster of close eigenvalues, select a new -* shift close to the cluster, find a new factorization, and refine -* the shifted eigenvalues to suitable accuracy. -* (d) For each eigenvalue with a large enough relative separation compute -* the corresponding eigenvector by forming a rank revealing twisted -* factorization. Go back to (c) for any clusters that remain. -* -* For more details, see: -* - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations -* to compute orthogonal eigenvectors of symmetric tridiagonal matrices," -* Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. -* - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and -* Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, -* 2004. Also LAPACK Working Note 154. -* - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric -* tridiagonal eigenvalue/eigenvector problem", -* Computer Science Division Technical Report No. UCB/CSD-97-971, -* UC Berkeley, May 1997. -* -* Further Details -* 1.ZSTEMR works only on machines which follow IEEE-754 -* floating-point standard in their handling of infinities and NaNs. -* This permits the use of efficient inner loops avoiding a check for -* zero divisors. -* -* 2. LAPACK routines can be used to reduce a complex Hermitean matrix to -* real symmetric tridiagonal form. -* -* (Any complex Hermitean tridiagonal matrix has real values on its diagonal -* and potentially complex numbers on its off-diagonals. By applying a -* similarity transform with an appropriate diagonal matrix -* diag(1,e^{i \phy_1}, ... , e^{i \phy_{n-1}}), the complex Hermitean -* matrix can be transformed into a real symmetric matrix and complex -* arithmetic can be entirely avoided.) -* -* While the eigenvectors of the real symmetric tridiagonal matrix are real, -* the eigenvectors of original complex Hermitean matrix have complex entries -* in general. -* Since LAPACK drivers overwrite the matrix data with the eigenvectors, -* ZSTEMR accepts complex workspace to facilitate interoperability -* with ZUNMTR or ZUPMTR. -* -* Arguments -* ========= -* -* JOBZ (input) CHARACTER*1 -* = 'N': Compute eigenvalues only; -* = 'V': Compute eigenvalues and eigenvectors. -* -* RANGE (input) CHARACTER*1 -* = 'A': all eigenvalues will be found. -* = 'V': all eigenvalues in the half-open interval (VL,VU] -* will be found. -* = 'I': the IL-th through IU-th eigenvalues will be found. -* -* N (input) INTEGER -* The order of the matrix. N >= 0. -* -* D (input/output) DOUBLE PRECISION array, dimension (N) -* On entry, the N diagonal elements of the tridiagonal matrix -* T. On exit, D is overwritten. -* -* E (input/output) DOUBLE PRECISION array, dimension (N) -* On entry, the (N-1) subdiagonal elements of the tridiagonal -* matrix T in elements 1 to N-1 of E. E(N) need not be set on -* input, but is used internally as workspace. -* On exit, E is overwritten. -* -* VL (input) DOUBLE PRECISION -* VU (input) DOUBLE PRECISION -* If RANGE='V', the lower and upper bounds of the interval to -* be searched for eigenvalues. VL < VU. -* Not referenced if RANGE = 'A' or 'I'. -* -* IL (input) INTEGER -* IU (input) INTEGER -* If RANGE='I', the indices (in ascending order) of the -* smallest and largest eigenvalues to be returned. -* 1 <= IL <= IU <= N, if N > 0. -* Not referenced if RANGE = 'A' or 'V'. -* -* M (output) INTEGER -* The total number of eigenvalues found. 0 <= M <= N. -* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. -* -* W (output) DOUBLE PRECISION array, dimension (N) -* The first M elements contain the selected eigenvalues in -* ascending order. -* -* Z (output) COMPLEX*16 array, dimension (LDZ, max(1,M) ) -* If JOBZ = 'V', and if INFO = 0, then the first M columns of Z -* contain the orthonormal eigenvectors of the matrix T -* corresponding to the selected eigenvalues, with the i-th -* column of Z holding the eigenvector associated with W(i). -* If JOBZ = 'N', then Z is not referenced. -* Note: the user must ensure that at least max(1,M) columns are -* supplied in the array Z; if RANGE = 'V', the exact value of M -* is not known in advance and can be computed with a workspace -* query by setting NZC = -1, see below. -* -* LDZ (input) INTEGER -* The leading dimension of the array Z. LDZ >= 1, and if -* JOBZ = 'V', then LDZ >= max(1,N). -* -* NZC (input) INTEGER -* The number of eigenvectors to be held in the array Z. -* If RANGE = 'A', then NZC >= max(1,N). -* If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU]. -* If RANGE = 'I', then NZC >= IU-IL+1. -* If NZC = -1, then a workspace query is assumed; the -* routine calculates the number of columns of the array Z that -* are needed to hold the eigenvectors. -* This value is returned as the first entry of the Z array, and -* no error message related to NZC is issued by XERBLA. -* -* ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) -* The support of the eigenvectors in Z, i.e., the indices -* indicating the nonzero elements in Z. The i-th computed eigenvector -* is nonzero only in elements ISUPPZ( 2*i-1 ) through -* ISUPPZ( 2*i ). This is relevant in the case when the matrix -* is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0. -* -* TRYRAC (input/output) LOGICAL -* If TRYRAC.EQ..TRUE., indicates that the code should check whether -* the tridiagonal matrix defines its eigenvalues to high relative -* accuracy. If so, the code uses relative-accuracy preserving -* algorithms that might be (a bit) slower depending on the matrix. -* If the matrix does not define its eigenvalues to high relative -* accuracy, the code can uses possibly faster algorithms. -* If TRYRAC.EQ..FALSE., the code is not required to guarantee -* relatively accurate eigenvalues and can use the fastest possible -* techniques. -* On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix -* does not define its eigenvalues to high relative accuracy. -* -* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) -* On exit, if INFO = 0, WORK(1) returns the optimal -* (and minimal) LWORK. -* -* LWORK (input) INTEGER -* The dimension of the array WORK. LWORK >= max(1,18*N) -* if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. -* If LWORK = -1, then a workspace query is assumed; the routine -* only calculates the optimal size of the WORK array, returns -* this value as the first entry of the WORK array, and no error -* message related to LWORK is issued by XERBLA. -* -* IWORK (workspace/output) INTEGER array, dimension (LIWORK) -* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. -* -* LIWORK (input) INTEGER -* The dimension of the array IWORK. LIWORK >= max(1,10*N) -* if the eigenvectors are desired, and LIWORK >= max(1,8*N) -* if only the eigenvalues are to be computed. -* If LIWORK = -1, then a workspace query is assumed; the -* routine only calculates the optimal size of the IWORK array, -* returns this value as the first entry of the IWORK array, and -* no error message related to LIWORK is issued by XERBLA. -* -* INFO (output) INTEGER -* On exit, INFO -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* > 0: if INFO = 1X, internal error in DLARRE, -* if INFO = 2X, internal error in ZLARRV. -* Here, the digit X = ABS( IINFO ) < 10, where IINFO is -* the nonzero error code returned by DLARRE or -* ZLARRV, respectively. -* -* -* Further Details -* =============== -* -* Based on contributions by -* Beresford Parlett, University of California, Berkeley, USA -* Jim Demmel, University of California, Berkeley, USA -* Inderjit Dhillon, University of Texas, Austin, USA -* Osni Marques, LBNL/NERSC, USA -* Christof Voemel, University of California, Berkeley, USA -* * ===================================================================== * * .. Parameters ..