zgeesx.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 409 行 · 第 1/3 页

</span><span class="comment">*</span><span class="comment">          that an error is only returned if LWORK &lt; max(1,2*N), but if
</span><span class="comment">*</span><span class="comment">          SENSE = 'E' or 'V' or 'B' this may not be large enough.
</span><span class="comment">*</span><span class="comment">          For good performance, LWORK must generally be larger.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates upper bound on the optimal size of the
</span><span class="comment">*</span><span class="comment">          array WORK, returns this value as the first entry of the WORK
</span><span class="comment">*</span><span class="comment">          array, and no error message related to LWORK is issued by
</span><span class="comment">*</span><span class="comment">          <a name="XERBLA.127"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  BWORK   (workspace) LOGICAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          Not referenced if SORT = 'N'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0: successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0: if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">          &gt; 0: if INFO = i, and i is
</span><span class="comment">*</span><span class="comment">             &lt;= N: the QR algorithm failed to compute all the
</span><span class="comment">*</span><span class="comment">                   eigenvalues; elements 1:ILO-1 and i+1:N of W
</span><span class="comment">*</span><span class="comment">                   contain those eigenvalues which have converged; if
</span><span class="comment">*</span><span class="comment">                   JOBVS = 'V', VS contains the transformation which
</span><span class="comment">*</span><span class="comment">                   reduces A to its partially converged Schur form.
</span><span class="comment">*</span><span class="comment">             = N+1: the eigenvalues could not be reordered because some
</span><span class="comment">*</span><span class="comment">                   eigenvalues were too close to separate (the problem
</span><span class="comment">*</span><span class="comment">                   is very ill-conditioned);
</span><span class="comment">*</span><span class="comment">             = N+2: after reordering, roundoff changed values of some
</span><span class="comment">*</span><span class="comment">                   complex eigenvalues so that leading eigenvalues in
</span><span class="comment">*</span><span class="comment">                   the Schur form no longer satisfy SELECT=.TRUE.  This
</span><span class="comment">*</span><span class="comment">                   could also be caused by underflow due to scaling.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            SCALEA, WANTSB, WANTSE, WANTSN, WANTST, WANTSV,
     $                   WANTVS
      INTEGER            HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
     $                   ITAU, IWRK, LWRK, MAXWRK, MINWRK
      DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SMLNUM
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      DOUBLE PRECISION   DUM( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLABAD.168"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, <a name="DLASCL.168"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="XERBLA.168"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, ZCOPY, <a name="ZGEBAK.168"></a><a href="zgebak.f.html#ZGEBAK.1">ZGEBAK</a>, <a name="ZGEBAL.168"></a><a href="zgebal.f.html#ZGEBAL.1">ZGEBAL</a>,
     $                   <a name="ZGEHRD.169"></a><a href="zgehrd.f.html#ZGEHRD.1">ZGEHRD</a>, <a name="ZHSEQR.169"></a><a href="zhseqr.f.html#ZHSEQR.1">ZHSEQR</a>, <a name="ZLACPY.169"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLASCL.169"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>, <a name="ZTRSEN.169"></a><a href="ztrsen.f.html#ZTRSEN.1">ZTRSEN</a>, <a name="ZUNGHR.169"></a><a href="zunghr.f.html#ZUNGHR.1">ZUNGHR</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.172"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.173"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      DOUBLE PRECISION   <a name="DLAMCH.174"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGE.174"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>
      EXTERNAL           <a name="LSAME.175"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.175"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="DLAMCH.175"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGE.175"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input arguments
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      WANTVS = <a name="LSAME.185"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVS, <span class="string">'V'</span> )
      WANTST = <a name="LSAME.186"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SORT, <span class="string">'S'</span> )
      WANTSN = <a name="LSAME.187"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SENSE, <span class="string">'N'</span> )
      WANTSE = <a name="LSAME.188"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SENSE, <span class="string">'E'</span> )
      WANTSV = <a name="LSAME.189"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SENSE, <span class="string">'V'</span> )
      WANTSB = <a name="LSAME.190"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SENSE, <span class="string">'B'</span> )
      IF( ( .NOT.WANTVS ) .AND. ( .NOT.<a name="LSAME.191"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVS, <span class="string">'N'</span> ) ) ) THEN
         INFO = -1
      ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.<a name="LSAME.193"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SORT, <span class="string">'N'</span> ) ) ) THEN
         INFO = -2
      ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
     $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
         INFO = -4
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
         INFO = -11
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute workspace
</span><span class="comment">*</span><span class="comment">      (Note: Comments in the code beginning &quot;Workspace:&quot; describe the
</span><span class="comment">*</span><span class="comment">       minimal amount of real workspace needed at that point in the
</span><span class="comment">*</span><span class="comment">       code, as well as the preferred amount for good performance.
</span><span class="comment">*</span><span class="comment">       CWorkspace refers to complex workspace, and RWorkspace to real
</span><span class="comment">*</span><span class="comment">       workspace. NB refers to the optimal block size for the
</span><span class="comment">*</span><span class="comment">       immediately following subroutine, as returned by <a name="ILAENV.212"></a><a href="hfy-index.html#ILAENV">ILAENV</a>.
</span><span class="comment">*</span><span class="comment">       HSWORK refers to the workspace preferred by <a name="ZHSEQR.213"></a><a href="zhseqr.f.html#ZHSEQR.1">ZHSEQR</a>, as
</span><span class="comment">*</span><span class="comment">       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
</span><span class="comment">*</span><span class="comment">       the worst case.
</span><span class="comment">*</span><span class="comment">       If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
</span><span class="comment">*</span><span class="comment">       depends on SDIM, which is computed by the routine <a name="ZTRSEN.217"></a><a href="ztrsen.f.html#ZTRSEN.1">ZTRSEN</a> later
</span><span class="comment">*</span><span class="comment">       in the code.)
</span><span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         IF( N.EQ.0 ) THEN
            MINWRK = 1
            LWRK = 1
         ELSE
            MAXWRK = N + N*<a name="ILAENV.225"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZGEHRD.225"></a><a href="zgehrd.f.html#ZGEHRD.1">ZGEHRD</a>'</span>, <span class="string">' '</span>, N, 1, N, 0 )
            MINWRK = 2*N
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="ZHSEQR.228"></a><a href="zhseqr.f.html#ZHSEQR.1">ZHSEQR</a>( <span class="string">'S'</span>, JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
     $             WORK, -1, IEVAL )
            HSWORK = WORK( 1 )
<span class="comment">*</span><span class="comment">
</span>            IF( .NOT.WANTVS ) THEN
               MAXWRK = MAX( MAXWRK, HSWORK )
            ELSE
               MAXWRK = MAX( MAXWRK, N + ( N - 1 )*<a name="ILAENV.235"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGHR.235"></a><a href="zunghr.f.html#ZUNGHR.1">ZUNGHR</a>'</span>,
     $                       <span class="string">' '</span>, N, 1, N, -1 ) )
               MAXWRK = MAX( MAXWRK, HSWORK )
            END IF
            LWRK = MAXWRK
            IF( .NOT.WANTSN )
     $         LWRK = MAX( LWRK, ( N*N )/2 )
         END IF
         WORK( 1 ) = LWRK
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.MINWRK ) THEN
            INFO = -15
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.251"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGEESX.251"></a><a href="zgeesx.f.html#ZGEESX.1">ZGEESX</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible