--- rpl/lapack/lapack/dtrsna.f 2010/08/13 21:04:00 1.6
+++ 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,