--- rpl/lapack/lapack/dtrsna.f 2010/08/06 15:32:36 1.4 +++ rpl/lapack/lapack/dtrsna.f 2023/08/07 08:39:13 1.19 @@ -1,13 +1,271 @@ +*> \brief \b DTRSNA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DTRSNA + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, +* LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, +* INFO ) +* +* .. Scalar Arguments .. +* CHARACTER HOWMNY, JOB +* INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N +* .. +* .. Array Arguments .. +* LOGICAL SELECT( * ) +* INTEGER IWORK( * ) +* DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ), +* $ VR( LDVR, * ), WORK( LDWORK, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DTRSNA estimates reciprocal condition numbers for specified +*> eigenvalues and/or right eigenvectors of a real upper +*> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q +*> orthogonal). +*> +*> T must be in Schur canonical form (as returned by DHSEQR), that is, +*> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each +*> 2-by-2 diagonal block has its diagonal elements equal and its +*> off-diagonal elements of opposite sign. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOB +*> \verbatim +*> JOB is CHARACTER*1 +*> Specifies whether condition numbers are required for +*> eigenvalues (S) or eigenvectors (SEP): +*> = 'E': for eigenvalues only (S); +*> = 'V': for eigenvectors only (SEP); +*> = 'B': for both eigenvalues and eigenvectors (S and SEP). +*> \endverbatim +*> +*> \param[in] HOWMNY +*> \verbatim +*> HOWMNY is CHARACTER*1 +*> = 'A': compute condition numbers for all eigenpairs; +*> = 'S': compute condition numbers for selected eigenpairs +*> specified by the array SELECT. +*> \endverbatim +*> +*> \param[in] SELECT +*> \verbatim +*> SELECT is LOGICAL array, dimension (N) +*> If HOWMNY = 'S', SELECT specifies the eigenpairs for which +*> condition numbers are required. To select condition numbers +*> for the eigenpair corresponding to a real eigenvalue w(j), +*> SELECT(j) must be set to .TRUE.. To select condition numbers +*> corresponding to a complex conjugate pair of eigenvalues w(j) +*> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be +*> set to .TRUE.. +*> If HOWMNY = 'A', SELECT is not referenced. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix T. N >= 0. +*> \endverbatim +*> +*> \param[in] T +*> \verbatim +*> T is DOUBLE PRECISION array, dimension (LDT,N) +*> The upper quasi-triangular matrix T, in Schur canonical form. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= max(1,N). +*> \endverbatim +*> +*> \param[in] VL +*> \verbatim +*> VL is DOUBLE PRECISION array, dimension (LDVL,M) +*> If JOB = 'E' or 'B', VL must contain left eigenvectors of T +*> (or of any Q*T*Q**T with Q orthogonal), corresponding to the +*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors +*> must be stored in consecutive columns of VL, as returned by +*> DHSEIN or DTREVC. +*> If JOB = 'V', VL is not referenced. +*> \endverbatim +*> +*> \param[in] LDVL +*> \verbatim +*> LDVL is INTEGER +*> The leading dimension of the array VL. +*> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N. +*> \endverbatim +*> +*> \param[in] VR +*> \verbatim +*> VR is DOUBLE PRECISION array, dimension (LDVR,M) +*> If JOB = 'E' or 'B', VR must contain right eigenvectors of T +*> (or of any Q*T*Q**T with Q orthogonal), corresponding to the +*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors +*> must be stored in consecutive columns of VR, as returned by +*> DHSEIN or DTREVC. +*> If JOB = 'V', VR is not referenced. +*> \endverbatim +*> +*> \param[in] LDVR +*> \verbatim +*> LDVR is INTEGER +*> The leading dimension of the array VR. +*> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N. +*> \endverbatim +*> +*> \param[out] S +*> \verbatim +*> S is DOUBLE PRECISION array, dimension (MM) +*> If JOB = 'E' or 'B', the reciprocal condition numbers of the +*> selected eigenvalues, stored in consecutive elements of the +*> array. For a complex conjugate pair of eigenvalues two +*> consecutive elements of S are set to the same value. Thus +*> S(j), SEP(j), and the j-th columns of VL and VR all +*> correspond to the same eigenpair (but not in general the +*> j-th eigenpair, unless all eigenpairs are selected). +*> If JOB = 'V', S is not referenced. +*> \endverbatim +*> +*> \param[out] SEP +*> \verbatim +*> SEP is DOUBLE PRECISION array, dimension (MM) +*> If JOB = 'V' or 'B', the estimated reciprocal condition +*> numbers of the selected eigenvectors, stored in consecutive +*> elements of the array. For a complex eigenvector two +*> consecutive elements of SEP are set to the same value. If +*> the eigenvalues cannot be reordered to compute SEP(j), SEP(j) +*> is set to 0; this can only occur when the true value would be +*> very small anyway. +*> If JOB = 'E', SEP is not referenced. +*> \endverbatim +*> +*> \param[in] MM +*> \verbatim +*> MM is INTEGER +*> The number of elements in the arrays S (if JOB = 'E' or 'B') +*> and/or SEP (if JOB = 'V' or 'B'). MM >= M. +*> \endverbatim +*> +*> \param[out] M +*> \verbatim +*> M is INTEGER +*> The number of elements of the arrays S and/or SEP actually +*> used to store the estimated condition numbers. +*> If HOWMNY = 'A', M is set to N. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (LDWORK,N+6) +*> If JOB = 'E', WORK is not referenced. +*> \endverbatim +*> +*> \param[in] LDWORK +*> \verbatim +*> LDWORK is INTEGER +*> The leading dimension of the array WORK. +*> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*(N-1)) +*> If JOB = 'E', IWORK 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 +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup doubleOTHERcomputational +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The reciprocal of the condition number of an eigenvalue lambda is +*> defined as +*> +*> S(lambda) = |v**T*u| / (norm(u)*norm(v)) +*> +*> where u and v are the right and left eigenvectors of T corresponding +*> to lambda; v**T denotes the transpose of v, and norm(u) +*> denotes the Euclidean norm. These reciprocal condition numbers always +*> lie between zero (very badly conditioned) and one (very well +*> conditioned). If n = 1, S(lambda) is defined to be 1. +*> +*> An approximate error bound for a computed eigenvalue W(i) is given by +*> +*> EPS * norm(T) / S(i) +*> +*> where EPS is the machine precision. +*> +*> The reciprocal of the condition number of the right eigenvector u +*> corresponding to lambda is defined as follows. Suppose +*> +*> T = ( lambda c ) +*> ( 0 T22 ) +*> +*> Then the reciprocal condition number is +*> +*> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I ) +*> +*> where sigma-min denotes the smallest singular value. We approximate +*> the smallest singular value by the reciprocal of an estimate of the +*> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is +*> defined to be abs(T(1,1)). +*> +*> An approximate error bound for a computed right eigenvector VR(i) +*> is given by +*> +*> EPS * norm(T) / SEP(i) +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, $ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, $ 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 -* -* Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH. * * .. Scalar Arguments .. CHARACTER HOWMNY, JOB @@ -20,160 +278,6 @@ $ VR( LDVR, * ), WORK( LDWORK, * ) * .. * -* Purpose -* ======= -* -* DTRSNA estimates reciprocal condition numbers for specified -* eigenvalues and/or right eigenvectors of a real upper -* quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q -* orthogonal). -* -* T must be in Schur canonical form (as returned by DHSEQR), that is, -* block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each -* 2-by-2 diagonal block has its diagonal elements equal and its -* off-diagonal elements of opposite sign. -* -* Arguments -* ========= -* -* JOB (input) CHARACTER*1 -* Specifies whether condition numbers are required for -* eigenvalues (S) or eigenvectors (SEP): -* = 'E': for eigenvalues only (S); -* = 'V': for eigenvectors only (SEP); -* = 'B': for both eigenvalues and eigenvectors (S and SEP). -* -* HOWMNY (input) CHARACTER*1 -* = 'A': compute condition numbers for all eigenpairs; -* = 'S': compute condition numbers for selected eigenpairs -* specified by the array SELECT. -* -* SELECT (input) LOGICAL array, dimension (N) -* If HOWMNY = 'S', SELECT specifies the eigenpairs for which -* condition numbers are required. To select condition numbers -* for the eigenpair corresponding to a real eigenvalue w(j), -* SELECT(j) must be set to .TRUE.. To select condition numbers -* corresponding to a complex conjugate pair of eigenvalues w(j) -* and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be -* set to .TRUE.. -* If HOWMNY = 'A', SELECT is not referenced. -* -* N (input) INTEGER -* The order of the matrix T. N >= 0. -* -* T (input) DOUBLE PRECISION array, dimension (LDT,N) -* The upper quasi-triangular matrix T, in Schur canonical form. -* -* LDT (input) INTEGER -* The leading dimension of the array T. LDT >= max(1,N). -* -* VL (input) DOUBLE PRECISION array, dimension (LDVL,M) -* If JOB = 'E' or 'B', VL must contain left eigenvectors of T -* (or of any Q*T*Q**T with Q orthogonal), corresponding to the -* eigenpairs specified by HOWMNY and SELECT. The eigenvectors -* must be stored in consecutive columns of VL, as returned by -* DHSEIN or DTREVC. -* If JOB = 'V', VL is not referenced. -* -* LDVL (input) INTEGER -* The leading dimension of the array VL. -* LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N. -* -* VR (input) DOUBLE PRECISION array, dimension (LDVR,M) -* If JOB = 'E' or 'B', VR must contain right eigenvectors of T -* (or of any Q*T*Q**T with Q orthogonal), corresponding to the -* eigenpairs specified by HOWMNY and SELECT. The eigenvectors -* must be stored in consecutive columns of VR, as returned by -* DHSEIN or DTREVC. -* If JOB = 'V', VR is not referenced. -* -* LDVR (input) INTEGER -* The leading dimension of the array VR. -* LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N. -* -* S (output) DOUBLE PRECISION array, dimension (MM) -* If JOB = 'E' or 'B', the reciprocal condition numbers of the -* selected eigenvalues, stored in consecutive elements of the -* array. For a complex conjugate pair of eigenvalues two -* consecutive elements of S are set to the same value. Thus -* S(j), SEP(j), and the j-th columns of VL and VR all -* correspond to the same eigenpair (but not in general the -* j-th eigenpair, unless all eigenpairs are selected). -* If JOB = 'V', S is not referenced. -* -* SEP (output) DOUBLE PRECISION array, dimension (MM) -* If JOB = 'V' or 'B', the estimated reciprocal condition -* numbers of the selected eigenvectors, stored in consecutive -* elements of the array. For a complex eigenvector two -* consecutive elements of SEP are set to the same value. If -* the eigenvalues cannot be reordered to compute SEP(j), SEP(j) -* is set to 0; this can only occur when the true value would be -* very small anyway. -* If JOB = 'E', SEP is not referenced. -* -* MM (input) INTEGER -* The number of elements in the arrays S (if JOB = 'E' or 'B') -* and/or SEP (if JOB = 'V' or 'B'). MM >= M. -* -* M (output) INTEGER -* The number of elements of the arrays S and/or SEP actually -* used to store the estimated condition numbers. -* If HOWMNY = 'A', M is set to N. -* -* WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,N+6) -* If JOB = 'E', WORK is not referenced. -* -* LDWORK (input) INTEGER -* The leading dimension of the array WORK. -* LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N. -* -* IWORK (workspace) INTEGER array, dimension (2*(N-1)) -* If JOB = 'E', IWORK is not referenced. -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* -* Further Details -* =============== -* -* The reciprocal of the condition number of an eigenvalue lambda is -* defined as -* -* S(lambda) = |v'*u| / (norm(u)*norm(v)) -* -* where u and v are the right and left eigenvectors of T corresponding -* to lambda; v' denotes the conjugate-transpose of v, and norm(u) -* denotes the Euclidean norm. These reciprocal condition numbers always -* lie between zero (very badly conditioned) and one (very well -* conditioned). If n = 1, S(lambda) is defined to be 1. -* -* An approximate error bound for a computed eigenvalue W(i) is given by -* -* EPS * norm(T) / S(i) -* -* where EPS is the machine precision. -* -* The reciprocal of the condition number of the right eigenvector u -* corresponding to lambda is defined as follows. Suppose -* -* T = ( lambda c ) -* ( 0 T22 ) -* -* Then the reciprocal condition number is -* -* SEP( lambda, T22 ) = sigma-min( T22 - lambda*I ) -* -* where sigma-min denotes the smallest singular value. We approximate -* the smallest singular value by the reciprocal of an estimate of the -* one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is -* defined to be abs(T(1,1)). -* -* An approximate error bound for a computed right eigenvector VR(i) -* is given by -* -* EPS * norm(T) / SEP(i) -* * ===================================================================== * * .. Parameters .. @@ -196,7 +300,7 @@ EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2 * .. * .. External Subroutines .. - EXTERNAL DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA + EXTERNAL DLABAD, DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT @@ -403,12 +507,12 @@ * * Form * -* C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ] -* [ mu ] -* [ .. ] -* [ .. ] -* [ mu ] -* where C' is conjugate transpose of complex matrix C, +* C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ] +* [ mu ] +* [ .. ] +* [ .. ] +* [ mu ] +* where C**T is transpose of matrix C, * and RWORK is stored starting in the N+1-st column of * WORK. * @@ -426,7 +530,7 @@ NN = 2*( N-1 ) END IF * -* Estimate norm(inv(C')) +* Estimate norm(inv(C**T)) * EST = ZERO KASE = 0 @@ -437,7 +541,7 @@ IF( KASE.EQ.1 ) THEN IF( N2.EQ.1 ) THEN * -* Real eigenvalue: solve C'*x = scale*c. +* Real eigenvalue: solve C**T*x = scale*c. * CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ), $ LDWORK, DUMMY, DUMM, SCALE, @@ -446,7 +550,7 @@ ELSE * * Complex eigenvalue: solve -* C'*(p+iq) = scale*(c+id) in real arithmetic. +* C**T*(p+iq) = scale*(c+id) in real arithmetic. * CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ), $ LDWORK, WORK( 1, N+1 ), MU, SCALE,