version 1.6, 2010/08/13 21:03:47
|
version 1.16, 2017/06/17 11:06:20
|
Line 1
|
Line 1
|
|
*> \brief \b DHSEIN |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DHSEIN + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhsein.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhsein.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhsein.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI, |
|
* VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL, |
|
* IFAILR, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER EIGSRC, INITV, SIDE |
|
* INTEGER INFO, LDH, LDVL, LDVR, M, MM, N |
|
* .. |
|
* .. Array Arguments .. |
|
* LOGICAL SELECT( * ) |
|
* INTEGER IFAILL( * ), IFAILR( * ) |
|
* DOUBLE PRECISION H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ), |
|
* $ WI( * ), WORK( * ), WR( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DHSEIN uses inverse iteration to find specified right and/or left |
|
*> eigenvectors of a real upper Hessenberg matrix H. |
|
*> |
|
*> The right eigenvector x and the left eigenvector y of the matrix H |
|
*> corresponding to an eigenvalue w are defined by: |
|
*> |
|
*> H * x = w * x, y**h * H = w * y**h |
|
*> |
|
*> where y**h denotes the conjugate transpose of the vector y. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] SIDE |
|
*> \verbatim |
|
*> SIDE is CHARACTER*1 |
|
*> = 'R': compute right eigenvectors only; |
|
*> = 'L': compute left eigenvectors only; |
|
*> = 'B': compute both right and left eigenvectors. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] EIGSRC |
|
*> \verbatim |
|
*> EIGSRC is CHARACTER*1 |
|
*> Specifies the source of eigenvalues supplied in (WR,WI): |
|
*> = 'Q': the eigenvalues were found using DHSEQR; thus, if |
|
*> H has zero subdiagonal elements, and so is |
|
*> block-triangular, then the j-th eigenvalue can be |
|
*> assumed to be an eigenvalue of the block containing |
|
*> the j-th row/column. This property allows DHSEIN to |
|
*> perform inverse iteration on just one diagonal block. |
|
*> = 'N': no assumptions are made on the correspondence |
|
*> between eigenvalues and diagonal blocks. In this |
|
*> case, DHSEIN must always perform inverse iteration |
|
*> using the whole matrix H. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] INITV |
|
*> \verbatim |
|
*> INITV is CHARACTER*1 |
|
*> = 'N': no initial vectors are supplied; |
|
*> = 'U': user-supplied initial vectors are stored in the arrays |
|
*> VL and/or VR. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] SELECT |
|
*> \verbatim |
|
*> SELECT is LOGICAL array, dimension (N) |
|
*> Specifies the eigenvectors to be computed. To select the |
|
*> real eigenvector corresponding to a real eigenvalue WR(j), |
|
*> SELECT(j) must be set to .TRUE.. To select the complex |
|
*> eigenvector corresponding to a complex eigenvalue |
|
*> (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)), |
|
*> either SELECT(j) or SELECT(j+1) or both must be set to |
|
*> .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is |
|
*> .FALSE.. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The order of the matrix H. N >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] H |
|
*> \verbatim |
|
*> H is DOUBLE PRECISION array, dimension (LDH,N) |
|
*> The upper Hessenberg matrix H. |
|
*> If a NaN is detected in H, the routine will return with INFO=-6. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDH |
|
*> \verbatim |
|
*> LDH is INTEGER |
|
*> The leading dimension of the array H. LDH >= max(1,N). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] WR |
|
*> \verbatim |
|
*> WR is DOUBLE PRECISION array, dimension (N) |
|
*> \endverbatim |
|
*> |
|
*> \param[in] WI |
|
*> \verbatim |
|
*> WI is DOUBLE PRECISION array, dimension (N) |
|
*> |
|
*> On entry, the real and imaginary parts of the eigenvalues of |
|
*> H; a complex conjugate pair of eigenvalues must be stored in |
|
*> consecutive elements of WR and WI. |
|
*> On exit, WR may have been altered since close eigenvalues |
|
*> are perturbed slightly in searching for independent |
|
*> eigenvectors. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] VL |
|
*> \verbatim |
|
*> VL is DOUBLE PRECISION array, dimension (LDVL,MM) |
|
*> On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must |
|
*> contain starting vectors for the inverse iteration for the |
|
*> left eigenvectors; the starting vector for each eigenvector |
|
*> must be in the same column(s) in which the eigenvector will |
|
*> be stored. |
|
*> On exit, if SIDE = 'L' or 'B', the left eigenvectors |
|
*> specified by SELECT will be stored consecutively in the |
|
*> columns of VL, in the same order as their eigenvalues. A |
|
*> complex eigenvector corresponding to a complex eigenvalue is |
|
*> stored in two consecutive columns, the first holding the real |
|
*> part and the second the imaginary part. |
|
*> If SIDE = 'R', VL is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDVL |
|
*> \verbatim |
|
*> LDVL is INTEGER |
|
*> The leading dimension of the array VL. |
|
*> LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] VR |
|
*> \verbatim |
|
*> VR is DOUBLE PRECISION array, dimension (LDVR,MM) |
|
*> On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must |
|
*> contain starting vectors for the inverse iteration for the |
|
*> right eigenvectors; the starting vector for each eigenvector |
|
*> must be in the same column(s) in which the eigenvector will |
|
*> be stored. |
|
*> On exit, if SIDE = 'R' or 'B', the right eigenvectors |
|
*> specified by SELECT will be stored consecutively in the |
|
*> columns of VR, in the same order as their eigenvalues. A |
|
*> complex eigenvector corresponding to a complex eigenvalue is |
|
*> stored in two consecutive columns, the first holding the real |
|
*> part and the second the imaginary part. |
|
*> If SIDE = 'L', VR is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDVR |
|
*> \verbatim |
|
*> LDVR is INTEGER |
|
*> The leading dimension of the array VR. |
|
*> LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] MM |
|
*> \verbatim |
|
*> MM is INTEGER |
|
*> The number of columns in the arrays VL and/or VR. MM >= M. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The number of columns in the arrays VL and/or VR required to |
|
*> store the eigenvectors; each selected real eigenvector |
|
*> occupies one column and each selected complex eigenvector |
|
*> occupies two columns. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension ((N+2)*N) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IFAILL |
|
*> \verbatim |
|
*> IFAILL is INTEGER array, dimension (MM) |
|
*> If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left |
|
*> eigenvector in the i-th column of VL (corresponding to the |
|
*> eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the |
|
*> eigenvector converged satisfactorily. If the i-th and (i+1)th |
|
*> columns of VL hold a complex eigenvector, then IFAILL(i) and |
|
*> IFAILL(i+1) are set to the same value. |
|
*> If SIDE = 'R', IFAILL is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IFAILR |
|
*> \verbatim |
|
*> IFAILR is INTEGER array, dimension (MM) |
|
*> If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right |
|
*> eigenvector in the i-th column of VR (corresponding to the |
|
*> eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the |
|
*> eigenvector converged satisfactorily. If the i-th and (i+1)th |
|
*> columns of VR hold a complex eigenvector, then IFAILR(i) and |
|
*> IFAILR(i+1) are set to the same value. |
|
*> If SIDE = 'L', IFAILR is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value |
|
*> > 0: if INFO = i, i is the number of eigenvectors which |
|
*> failed to converge; see IFAILL and IFAILR for further |
|
*> details. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date December 2016 |
|
* |
|
*> \ingroup doubleOTHERcomputational |
|
* |
|
*> \par Further Details: |
|
* ===================== |
|
*> |
|
*> \verbatim |
|
*> |
|
*> Each eigenvector is normalized so that the element of largest |
|
*> magnitude has magnitude 1; here the magnitude of a complex number |
|
*> (x,y) is taken to be |x|+|y|. |
|
*> \endverbatim |
|
*> |
|
* ===================================================================== |
SUBROUTINE DHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI, |
SUBROUTINE DHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI, |
$ VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL, |
$ VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL, |
$ IFAILR, INFO ) |
$ IFAILR, INFO ) |
* |
* |
* -- LAPACK routine (version 3.2) -- |
* -- LAPACK computational routine (version 3.7.0) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2006 |
* December 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER EIGSRC, INITV, SIDE |
CHARACTER EIGSRC, INITV, SIDE |
Line 18
|
Line 279
|
$ WI( * ), WORK( * ), WR( * ) |
$ WI( * ), WORK( * ), WR( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DHSEIN uses inverse iteration to find specified right and/or left |
|
* eigenvectors of a real upper Hessenberg matrix H. |
|
* |
|
* The right eigenvector x and the left eigenvector y of the matrix H |
|
* corresponding to an eigenvalue w are defined by: |
|
* |
|
* H * x = w * x, y**h * H = w * y**h |
|
* |
|
* where y**h denotes the conjugate transpose of the vector y. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* SIDE (input) CHARACTER*1 |
|
* = 'R': compute right eigenvectors only; |
|
* = 'L': compute left eigenvectors only; |
|
* = 'B': compute both right and left eigenvectors. |
|
* |
|
* EIGSRC (input) CHARACTER*1 |
|
* Specifies the source of eigenvalues supplied in (WR,WI): |
|
* = 'Q': the eigenvalues were found using DHSEQR; thus, if |
|
* H has zero subdiagonal elements, and so is |
|
* block-triangular, then the j-th eigenvalue can be |
|
* assumed to be an eigenvalue of the block containing |
|
* the j-th row/column. This property allows DHSEIN to |
|
* perform inverse iteration on just one diagonal block. |
|
* = 'N': no assumptions are made on the correspondence |
|
* between eigenvalues and diagonal blocks. In this |
|
* case, DHSEIN must always perform inverse iteration |
|
* using the whole matrix H. |
|
* |
|
* INITV (input) CHARACTER*1 |
|
* = 'N': no initial vectors are supplied; |
|
* = 'U': user-supplied initial vectors are stored in the arrays |
|
* VL and/or VR. |
|
* |
|
* SELECT (input/output) LOGICAL array, dimension (N) |
|
* Specifies the eigenvectors to be computed. To select the |
|
* real eigenvector corresponding to a real eigenvalue WR(j), |
|
* SELECT(j) must be set to .TRUE.. To select the complex |
|
* eigenvector corresponding to a complex eigenvalue |
|
* (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)), |
|
* either SELECT(j) or SELECT(j+1) or both must be set to |
|
* .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is |
|
* .FALSE.. |
|
* |
|
* N (input) INTEGER |
|
* The order of the matrix H. N >= 0. |
|
* |
|
* H (input) DOUBLE PRECISION array, dimension (LDH,N) |
|
* The upper Hessenberg matrix H. |
|
* |
|
* LDH (input) INTEGER |
|
* The leading dimension of the array H. LDH >= max(1,N). |
|
* |
|
* WR (input/output) DOUBLE PRECISION array, dimension (N) |
|
* WI (input) DOUBLE PRECISION array, dimension (N) |
|
* On entry, the real and imaginary parts of the eigenvalues of |
|
* H; a complex conjugate pair of eigenvalues must be stored in |
|
* consecutive elements of WR and WI. |
|
* On exit, WR may have been altered since close eigenvalues |
|
* are perturbed slightly in searching for independent |
|
* eigenvectors. |
|
* |
|
* VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) |
|
* On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must |
|
* contain starting vectors for the inverse iteration for the |
|
* left eigenvectors; the starting vector for each eigenvector |
|
* must be in the same column(s) in which the eigenvector will |
|
* be stored. |
|
* On exit, if SIDE = 'L' or 'B', the left eigenvectors |
|
* specified by SELECT will be stored consecutively in the |
|
* columns of VL, in the same order as their eigenvalues. A |
|
* complex eigenvector corresponding to a complex eigenvalue is |
|
* stored in two consecutive columns, the first holding the real |
|
* part and the second the imaginary part. |
|
* If SIDE = 'R', VL is not referenced. |
|
* |
|
* LDVL (input) INTEGER |
|
* The leading dimension of the array VL. |
|
* LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. |
|
* |
|
* VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) |
|
* On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must |
|
* contain starting vectors for the inverse iteration for the |
|
* right eigenvectors; the starting vector for each eigenvector |
|
* must be in the same column(s) in which the eigenvector will |
|
* be stored. |
|
* On exit, if SIDE = 'R' or 'B', the right eigenvectors |
|
* specified by SELECT will be stored consecutively in the |
|
* columns of VR, in the same order as their eigenvalues. A |
|
* complex eigenvector corresponding to a complex eigenvalue is |
|
* stored in two consecutive columns, the first holding the real |
|
* part and the second the imaginary part. |
|
* If SIDE = 'L', VR is not referenced. |
|
* |
|
* LDVR (input) INTEGER |
|
* The leading dimension of the array VR. |
|
* LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. |
|
* |
|
* MM (input) INTEGER |
|
* The number of columns in the arrays VL and/or VR. MM >= M. |
|
* |
|
* M (output) INTEGER |
|
* The number of columns in the arrays VL and/or VR required to |
|
* store the eigenvectors; each selected real eigenvector |
|
* occupies one column and each selected complex eigenvector |
|
* occupies two columns. |
|
* |
|
* WORK (workspace) DOUBLE PRECISION array, dimension ((N+2)*N) |
|
* |
|
* IFAILL (output) INTEGER array, dimension (MM) |
|
* If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left |
|
* eigenvector in the i-th column of VL (corresponding to the |
|
* eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the |
|
* eigenvector converged satisfactorily. If the i-th and (i+1)th |
|
* columns of VL hold a complex eigenvector, then IFAILL(i) and |
|
* IFAILL(i+1) are set to the same value. |
|
* If SIDE = 'R', IFAILL is not referenced. |
|
* |
|
* IFAILR (output) INTEGER array, dimension (MM) |
|
* If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right |
|
* eigenvector in the i-th column of VR (corresponding to the |
|
* eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the |
|
* eigenvector converged satisfactorily. If the i-th and (i+1)th |
|
* columns of VR hold a complex eigenvector, then IFAILR(i) and |
|
* IFAILR(i+1) are set to the same value. |
|
* If SIDE = 'L', IFAILR is not referenced. |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit |
|
* < 0: if INFO = -i, the i-th argument had an illegal value |
|
* > 0: if INFO = i, i is the number of eigenvectors which |
|
* failed to converge; see IFAILL and IFAILR for further |
|
* details. |
|
* |
|
* Further Details |
|
* =============== |
|
* |
|
* Each eigenvector is normalized so that the element of largest |
|
* magnitude has magnitude 1; here the magnitude of a complex number |
|
* (x,y) is taken to be |x|+|y|. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 177
|
Line 292
|
$ WKR |
$ WKR |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
LOGICAL LSAME |
LOGICAL LSAME, DISNAN |
DOUBLE PRECISION DLAMCH, DLANHS |
DOUBLE PRECISION DLAMCH, DLANHS |
EXTERNAL LSAME, DLAMCH, DLANHS |
EXTERNAL LSAME, DLAMCH, DLANHS, DISNAN |
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DLAEIN, XERBLA |
EXTERNAL DLAEIN, XERBLA |
Line 309
|
Line 424
|
* has not ben computed before. |
* has not ben computed before. |
* |
* |
HNORM = DLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, WORK ) |
HNORM = DLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, WORK ) |
IF( HNORM.GT.ZERO ) THEN |
IF( DISNAN( HNORM ) ) THEN |
|
INFO = -6 |
|
RETURN |
|
ELSE IF( HNORM.GT.ZERO ) THEN |
EPS3 = HNORM*ULP |
EPS3 = HNORM*ULP |
ELSE |
ELSE |
EPS3 = SMLNUM |
EPS3 = SMLNUM |