Diff for /rpl/lapack/lapack/zgeev.f between versions 1.5 and 1.18

version 1.5, 2010/08/07 13:22:30 version 1.18, 2020/05/21 21:46:03
Line 1 Line 1
   *> \brief <b> ZGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download ZGEEV + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeev.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeev.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeev.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
   *                         WORK, LWORK, RWORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          JOBVL, JOBVR
   *       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   RWORK( * )
   *       COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
   *      $                   W( * ), WORK( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
   *> eigenvalues and, optionally, the left and/or right eigenvectors.
   *>
   *> The right eigenvector v(j) of A satisfies
   *>                  A * v(j) = lambda(j) * v(j)
   *> where lambda(j) is its eigenvalue.
   *> The left eigenvector u(j) of A satisfies
   *>               u(j)**H * A = lambda(j) * u(j)**H
   *> where u(j)**H denotes the conjugate transpose of u(j).
   *>
   *> The computed eigenvectors are normalized to have Euclidean norm
   *> equal to 1 and largest component real.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] JOBVL
   *> \verbatim
   *>          JOBVL is CHARACTER*1
   *>          = 'N': left eigenvectors of A are not computed;
   *>          = 'V': left eigenvectors of are computed.
   *> \endverbatim
   *>
   *> \param[in] JOBVR
   *> \verbatim
   *>          JOBVR is CHARACTER*1
   *>          = 'N': right eigenvectors of A are not computed;
   *>          = 'V': right eigenvectors of A are computed.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the matrix A. N >= 0.
   *> \endverbatim
   *>
   *> \param[in,out] A
   *> \verbatim
   *>          A is COMPLEX*16 array, dimension (LDA,N)
   *>          On entry, the N-by-N matrix A.
   *>          On exit, A has been overwritten.
   *> \endverbatim
   *>
   *> \param[in] LDA
   *> \verbatim
   *>          LDA is INTEGER
   *>          The leading dimension of the array A.  LDA >= max(1,N).
   *> \endverbatim
   *>
   *> \param[out] W
   *> \verbatim
   *>          W is COMPLEX*16 array, dimension (N)
   *>          W contains the computed eigenvalues.
   *> \endverbatim
   *>
   *> \param[out] VL
   *> \verbatim
   *>          VL is COMPLEX*16 array, dimension (LDVL,N)
   *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
   *>          after another in the columns of VL, in the same order
   *>          as their eigenvalues.
   *>          If JOBVL = 'N', VL is not referenced.
   *>          u(j) = VL(:,j), the j-th column of VL.
   *> \endverbatim
   *>
   *> \param[in] LDVL
   *> \verbatim
   *>          LDVL is INTEGER
   *>          The leading dimension of the array VL.  LDVL >= 1; if
   *>          JOBVL = 'V', LDVL >= N.
   *> \endverbatim
   *>
   *> \param[out] VR
   *> \verbatim
   *>          VR is COMPLEX*16 array, dimension (LDVR,N)
   *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
   *>          after another in the columns of VR, in the same order
   *>          as their eigenvalues.
   *>          If JOBVR = 'N', VR is not referenced.
   *>          v(j) = VR(:,j), the j-th column of VR.
   *> \endverbatim
   *>
   *> \param[in] LDVR
   *> \verbatim
   *>          LDVR is INTEGER
   *>          The leading dimension of the array VR.  LDVR >= 1; if
   *>          JOBVR = 'V', LDVR >= N.
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
   *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   *> \endverbatim
   *>
   *> \param[in] LWORK
   *> \verbatim
   *>          LWORK is INTEGER
   *>          The dimension of the array WORK.  LWORK >= max(1,2*N).
   *>          For good performance, LWORK must generally be larger.
   *>
   *>          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] RWORK
   *> \verbatim
   *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
   *> \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, the QR algorithm failed to compute all the
   *>                eigenvalues, and no eigenvectors have been computed;
   *>                elements i+1:N of W contain eigenvalues which have
   *>                converged.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \date June 2016
   *
   *  @precisions fortran z -> c
   *
   *> \ingroup complex16GEeigen
   *
   *  =====================================================================
       SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,        SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
      $                  WORK, LWORK, RWORK, INFO )       $                  WORK, LWORK, RWORK, INFO )
         implicit none
 *  *
 *  -- LAPACK driver routine (version 3.2) --  *  -- LAPACK driver 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  *     June 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBVL, JOBVR        CHARACTER          JOBVL, JOBVR
Line 16 Line 195
      $                   W( * ), WORK( * )       $                   W( * ), WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the  
 *  eigenvalues and, optionally, the left and/or right eigenvectors.  
 *  
 *  The right eigenvector v(j) of A satisfies  
 *                   A * v(j) = lambda(j) * v(j)  
 *  where lambda(j) is its eigenvalue.  
 *  The left eigenvector u(j) of A satisfies  
 *                u(j)**H * A = lambda(j) * u(j)**H  
 *  where u(j)**H denotes the conjugate transpose of u(j).  
 *  
 *  The computed eigenvectors are normalized to have Euclidean norm  
 *  equal to 1 and largest component real.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  JOBVL   (input) CHARACTER*1  
 *          = 'N': left eigenvectors of A are not computed;  
 *          = 'V': left eigenvectors of are computed.  
 *  
 *  JOBVR   (input) CHARACTER*1  
 *          = 'N': right eigenvectors of A are not computed;  
 *          = 'V': right eigenvectors of A are computed.  
 *  
 *  N       (input) INTEGER  
 *          The order of the matrix A. N >= 0.  
 *  
 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)  
 *          On entry, the N-by-N matrix A.  
 *          On exit, A has been overwritten.  
 *  
 *  LDA     (input) INTEGER  
 *          The leading dimension of the array A.  LDA >= max(1,N).  
 *  
 *  W       (output) COMPLEX*16 array, dimension (N)  
 *          W contains the computed eigenvalues.  
 *  
 *  VL      (output) COMPLEX*16 array, dimension (LDVL,N)  
 *          If JOBVL = 'V', the left eigenvectors u(j) are stored one  
 *          after another in the columns of VL, in the same order  
 *          as their eigenvalues.  
 *          If JOBVL = 'N', VL is not referenced.  
 *          u(j) = VL(:,j), the j-th column of VL.  
 *  
 *  LDVL    (input) INTEGER  
 *          The leading dimension of the array VL.  LDVL >= 1; if  
 *          JOBVL = 'V', LDVL >= N.  
 *  
 *  VR      (output) COMPLEX*16 array, dimension (LDVR,N)  
 *          If JOBVR = 'V', the right eigenvectors v(j) are stored one  
 *          after another in the columns of VR, in the same order  
 *          as their eigenvalues.  
 *          If JOBVR = 'N', VR is not referenced.  
 *          v(j) = VR(:,j), the j-th column of VR.  
 *  
 *  LDVR    (input) INTEGER  
 *          The leading dimension of the array VR.  LDVR >= 1; if  
 *          JOBVR = 'V', LDVR >= N.  
 *  
 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))  
 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.  
 *  
 *  LWORK   (input) INTEGER  
 *          The dimension of the array WORK.  LWORK >= max(1,2*N).  
 *          For good performance, LWORK must generally be larger.  
 *  
 *          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.  
 *  
 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit  
 *          < 0:  if INFO = -i, the i-th argument had an illegal value.  
 *          > 0:  if INFO = i, the QR algorithm failed to compute all the  
 *                eigenvalues, and no eigenvectors have been computed;  
 *                elements and i+1:N of W contain eigenvalues which have  
 *                converged.  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 110 Line 205
       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR        LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
       CHARACTER          SIDE        CHARACTER          SIDE
       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,        INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
      $                   IWRK, K, MAXWRK, MINWRK, NOUT       $                   IWRK, K, LWORK_TREVC, MAXWRK, MINWRK, NOUT
       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM        DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
       COMPLEX*16         TMP        COMPLEX*16         TMP
 *     ..  *     ..
Line 120 Line 215
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           DLABAD, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD,        EXTERNAL           DLABAD, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD,
      $                   ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC, ZUNGHR       $                   ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC3, ZUNGHR
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME
Line 129 Line 224
       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE        EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT        INTRINSIC          DBLE, DCMPLX, CONJG, AIMAG, MAX, SQRT
 *     ..  *     ..
 *     .. Executable Statements ..  *     .. Executable Statements ..
 *  *
Line 174 Line 269
             IF( WANTVL ) THEN              IF( WANTVL ) THEN
                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',                 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
      $                       ' ', N, 1, N, -1 ) )       $                       ' ', N, 1, N, -1 ) )
                  CALL ZTREVC3( 'L', 'B', SELECT, N, A, LDA,
        $                       VL, LDVL, VR, LDVR,
        $                       N, NOUT, WORK, -1, RWORK, -1, IERR )
                  LWORK_TREVC = INT( WORK(1) )
                  MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,                 CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
      $                WORK, -1, INFO )       $                      WORK, -1, INFO )
             ELSE IF( WANTVR ) THEN              ELSE IF( WANTVR ) THEN
                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',                 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
      $                       ' ', N, 1, N, -1 ) )       $                       ' ', N, 1, N, -1 ) )
                  CALL ZTREVC3( 'R', 'B', SELECT, N, A, LDA,
        $                       VL, LDVL, VR, LDVR,
        $                       N, NOUT, WORK, -1, RWORK, -1, IERR )
                  LWORK_TREVC = INT( WORK(1) )
                  MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,                 CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
      $                WORK, -1, INFO )       $                      WORK, -1, INFO )
             ELSE              ELSE
                CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,                 CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
      $                WORK, -1, INFO )       $                      WORK, -1, INFO )
             END IF              END IF
             HSWORK = WORK( 1 )              HSWORK = INT( WORK(1) )
             MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )              MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
          END IF           END IF
          WORK( 1 ) = MAXWRK           WORK( 1 ) = MAXWRK
Line 312 Line 417
      $                WORK( IWRK ), LWORK-IWRK+1, INFO )       $                WORK( IWRK ), LWORK-IWRK+1, INFO )
       END IF        END IF
 *  *
 *     If INFO > 0 from ZHSEQR, then quit  *     If INFO .NE. 0 from ZHSEQR, then quit
 *  *
       IF( INFO.GT.0 )        IF( INFO.NE.0 )
      $   GO TO 50       $   GO TO 50
 *  *
       IF( WANTVL .OR. WANTVR ) THEN        IF( WANTVL .OR. WANTVR ) THEN
 *  *
 *        Compute left and/or right eigenvectors  *        Compute left and/or right eigenvectors
 *        (CWorkspace: need 2*N)  *        (CWorkspace: need 2*N, prefer N + 2*N*NB)
 *        (RWorkspace: need 2*N)  *        (RWorkspace: need 2*N)
 *  *
          IRWORK = IBAL + N           IRWORK = IBAL + N
          CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,           CALL ZTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
      $                N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )       $                 N, NOUT, WORK( IWRK ), LWORK-IWRK+1,
        $                 RWORK( IRWORK ), N, IERR )
       END IF        END IF
 *  *
       IF( WANTVL ) THEN        IF( WANTVL ) THEN
Line 344 Line 450
             CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )              CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
             DO 10 K = 1, N              DO 10 K = 1, N
                RWORK( IRWORK+K-1 ) = DBLE( VL( K, I ) )**2 +                 RWORK( IRWORK+K-1 ) = DBLE( VL( K, I ) )**2 +
      $                               DIMAG( VL( K, I ) )**2       $                               AIMAG( VL( K, I ) )**2
    10       CONTINUE     10       CONTINUE
             K = IDAMAX( N, RWORK( IRWORK ), 1 )              K = IDAMAX( N, RWORK( IRWORK ), 1 )
             TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )              TMP = CONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
             CALL ZSCAL( N, TMP, VL( 1, I ), 1 )              CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
             VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )              VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
    20    CONTINUE     20    CONTINUE
Line 369 Line 475
             CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )              CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
             DO 30 K = 1, N              DO 30 K = 1, N
                RWORK( IRWORK+K-1 ) = DBLE( VR( K, I ) )**2 +                 RWORK( IRWORK+K-1 ) = DBLE( VR( K, I ) )**2 +
      $                               DIMAG( VR( K, I ) )**2       $                               AIMAG( VR( K, I ) )**2
    30       CONTINUE     30       CONTINUE
             K = IDAMAX( N, RWORK( IRWORK ), 1 )              K = IDAMAX( N, RWORK( IRWORK ), 1 )
             TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )              TMP = CONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
             CALL ZSCAL( N, TMP, VR( 1, I ), 1 )              CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
             VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )              VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
    40    CONTINUE     40    CONTINUE

Removed from v.1.5  
changed lines
  Added in v.1.18


CVSweb interface <joel.bertrand@systella.fr>