[BACK]Return to dtrsna.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dtrsna.f between versions 1.4 and 1.19

version 1.4, 2010/08/06 15:32:36 version 1.19, 2023/08/07 08:39:13
Line 1 Line 1
   *> \brief \b DTRSNA
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DTRSNA + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna.f">
   *> [TXT]</a>
   *> \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,        SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,       $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
      $                   INFO )       $                   INFO )
 *  *
 *  -- LAPACK routine (version 3.2) --  *  -- LAPACK computational routine --
 *  -- 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  
 *  
 *     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          HOWMNY, JOB        CHARACTER          HOWMNY, JOB
Line 20 Line 278
      $                   VR( LDVR, * ), WORK( LDWORK, * )       $                   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 ..  *     .. Parameters ..
Line 196 Line 300
       EXTERNAL           LSAME, DDOT, DLAMCH, DLAPY2, DNRM2        EXTERNAL           LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA        EXTERNAL           DLABAD, DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          ABS, MAX, SQRT        INTRINSIC          ABS, MAX, SQRT
Line 403 Line 507
 *  *
 *                 Form  *                 Form
 *  *
 *                 C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]  *                 C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
 *                                        [   mu                     ]  *                                          [   mu                     ]
 *                                        [         ..               ]  *                                          [         ..               ]
 *                                        [             ..           ]  *                                          [             ..           ]
 *                                        [                  mu      ]  *                                          [                  mu      ]
 *                 where C' is conjugate transpose of complex matrix C,  *                 where C**T is transpose of matrix C,
 *                 and RWORK is stored starting in the N+1-st column of  *                 and RWORK is stored starting in the N+1-st column of
 *                 WORK.  *                 WORK.
 *  *
Line 426 Line 530
                   NN = 2*( N-1 )                    NN = 2*( N-1 )
                END IF                 END IF
 *  *
 *              Estimate norm(inv(C'))  *              Estimate norm(inv(C**T))
 *  *
                EST = ZERO                 EST = ZERO
                KASE = 0                 KASE = 0
Line 437 Line 541
                   IF( KASE.EQ.1 ) THEN                    IF( KASE.EQ.1 ) THEN
                      IF( N2.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 ),                          CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
      $                               LDWORK, DUMMY, DUMM, SCALE,       $                               LDWORK, DUMMY, DUMM, SCALE,
Line 446 Line 550
                      ELSE                       ELSE
 *  *
 *                       Complex eigenvalue: solve  *                       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 ),                          CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
      $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,       $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
Removed from v.1.4  
changed lines
  Added in v.1.19

CVSweb interface <joel.bertrand@systella.fr>