cgeesx.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) REAL 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>      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
<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
      REAL               ANRM, BIGNUM, CSCALE, EPS, SMLNUM
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      REAL               DUM( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CCOPY, <a name="CGEBAK.168"></a><a href="cgebak.f.html#CGEBAK.1">CGEBAK</a>, <a name="CGEBAL.168"></a><a href="cgebal.f.html#CGEBAL.1">CGEBAL</a>, <a name="CGEHRD.168"></a><a href="cgehrd.f.html#CGEHRD.1">CGEHRD</a>, <a name="CHSEQR.168"></a><a href="chseqr.f.html#CHSEQR.1">CHSEQR</a>, <a name="CLACPY.168"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>,
     $                   <a name="CLASCL.169"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>, <a name="CTRSEN.169"></a><a href="ctrsen.f.html#CTRSEN.1">CTRSEN</a>, <a name="CUNGHR.169"></a><a href="cunghr.f.html#CUNGHR.1">CUNGHR</a>, <a name="SLABAD.169"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>, <a name="SLASCL.169"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>, <a name="XERBLA.169"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</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>
      REAL               <a name="CLANGE.174"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>, <a name="SLAMCH.174"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</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="CLANGE.175"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>, <a name="SLAMCH.175"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</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="CHSEQR.213"></a><a href="chseqr.f.html#CHSEQR.1">CHSEQR</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="CTRSEN.217"></a><a href="ctrsen.f.html#CTRSEN.1">CTRSEN</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="CGEHRD.225"></a><a href="cgehrd.f.html#CGEHRD.1">CGEHRD</a>'</span>, <span class="string">' '</span>, N, 1, N, 0 )
            MINWRK = 2*N
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="CHSEQR.228"></a><a href="chseqr.f.html#CHSEQR.1">CHSEQR</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="CUNGHR.235"></a><a href="cunghr.f.html#CUNGHR.1">CUNGHR</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="CGEESX.251"></a><a href="cgeesx.f.html#CGEESX.1">CGEESX</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