dgegv.f.html

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

</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK.  LWORK &gt;= max(1,8*N).
</span><span class="comment">*</span><span class="comment">          For good performance, LWORK must generally be larger.
</span><span class="comment">*</span><span class="comment">          To compute the optimal value of LWORK, call <a name="ILAENV.159"></a><a href="hfy-index.html#ILAENV">ILAENV</a> to get
</span><span class="comment">*</span><span class="comment">          blocksizes (for <a name="DGEQRF.160"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>, <a name="DORMQR.160"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>, and <a name="DORGQR.160"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>.)  Then compute:
</span><span class="comment">*</span><span class="comment">          NB  -- MAX of the blocksizes for <a name="DGEQRF.161"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>, <a name="DORMQR.161"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>, and <a name="DORGQR.161"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>;
</span><span class="comment">*</span><span class="comment">          The optimal LWORK is:
</span><span class="comment">*</span><span class="comment">              2*N + MAX( 6*N, N*(NB+1) ).
</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 the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.168"></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">  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">          = 1,...,N:
</span><span class="comment">*</span><span class="comment">                The QZ iteration failed.  No eigenvectors have been
</span><span class="comment">*</span><span class="comment">                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
</span><span class="comment">*</span><span class="comment">                should be correct for j=INFO+1,...,N.
</span><span class="comment">*</span><span class="comment">          &gt; N:  errors that usually indicate LAPACK problems:
</span><span class="comment">*</span><span class="comment">                =N+1: error return from <a name="DGGBAL.178"></a><a href="dggbal.f.html#DGGBAL.1">DGGBAL</a>
</span><span class="comment">*</span><span class="comment">                =N+2: error return from <a name="DGEQRF.179"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>
</span><span class="comment">*</span><span class="comment">                =N+3: error return from <a name="DORMQR.180"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>
</span><span class="comment">*</span><span class="comment">                =N+4: error return from <a name="DORGQR.181"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>
</span><span class="comment">*</span><span class="comment">                =N+5: error return from <a name="DGGHRD.182"></a><a href="dgghrd.f.html#DGGHRD.1">DGGHRD</a>
</span><span class="comment">*</span><span class="comment">                =N+6: error return from <a name="DHGEQZ.183"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a> (other than failed
</span><span class="comment">*</span><span class="comment">                                                iteration)
</span><span class="comment">*</span><span class="comment">                =N+7: error return from <a name="DTGEVC.185"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>
</span><span class="comment">*</span><span class="comment">                =N+8: error return from <a name="DGGBAK.186"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a> (computing VL)
</span><span class="comment">*</span><span class="comment">                =N+9: error return from <a name="DGGBAK.187"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a> (computing VR)
</span><span class="comment">*</span><span class="comment">                =N+10: error return from <a name="DLASCL.188"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a> (various calls)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Balancing
</span><span class="comment">*</span><span class="comment">  ---------
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  This driver calls <a name="DGGBAL.196"></a><a href="dggbal.f.html#DGGBAL.1">DGGBAL</a> to both permute and scale rows and columns
</span><span class="comment">*</span><span class="comment">  of A and B.  The permutations PL and PR are chosen so that PL*A*PR
</span><span class="comment">*</span><span class="comment">  and PL*B*R will be upper triangular except for the diagonal blocks
</span><span class="comment">*</span><span class="comment">  A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
</span><span class="comment">*</span><span class="comment">  possible.  The diagonal scaling matrices DL and DR are chosen so
</span><span class="comment">*</span><span class="comment">  that the pair  DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
</span><span class="comment">*</span><span class="comment">  one (except for the elements that start out zero.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  After the eigenvalues and eigenvectors of the balanced matrices
</span><span class="comment">*</span><span class="comment">  have been computed, <a name="DGGBAK.205"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a> transforms the eigenvectors back to what
</span><span class="comment">*</span><span class="comment">  they would have been (in perfect arithmetic) if they had not been
</span><span class="comment">*</span><span class="comment">  balanced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Contents of A and B on Exit
</span><span class="comment">*</span><span class="comment">  -------- -- - --- - -- ----
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
</span><span class="comment">*</span><span class="comment">  both), then on exit the arrays A and B will contain the real Schur
</span><span class="comment">*</span><span class="comment">  form[*] of the &quot;balanced&quot; versions of A and B.  If no eigenvectors
</span><span class="comment">*</span><span class="comment">  are computed, then only the diagonal blocks will be correct.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  [*] See <a name="DHGEQZ.217"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>, <a name="DGEGS.217"></a><a href="dgegs.f.html#DGEGS.1">DGEGS</a>, or read the book &quot;Matrix Computations&quot;,
</span><span class="comment">*</span><span class="comment">      by Golub &amp; van Loan, pub. by Johns Hopkins U. Press.
</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            ILIMIT, ILV, ILVL, ILVR, LQUERY
      CHARACTER          CHTEMP
      INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
     $                   IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
     $                   LWKMIN, LWKOPT, NB, NB1, NB2, NB3
      DOUBLE PRECISION   ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
     $                   BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
     $                   SALFAI, SALFAR, SBETA, SCALE, TEMP
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      LOGICAL            LDUMMA( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DGEQRF.240"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>, <a name="DGGBAK.240"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a>, <a name="DGGBAL.240"></a><a href="dggbal.f.html#DGGBAL.1">DGGBAL</a>, <a name="DGGHRD.240"></a><a href="dgghrd.f.html#DGGHRD.1">DGGHRD</a>, <a name="DHGEQZ.240"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>, <a name="DLACPY.240"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>,
     $                   <a name="DLASCL.241"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="DLASET.241"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>, <a name="DORGQR.241"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>, <a name="DORMQR.241"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>, <a name="DTGEVC.241"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>, <a name="XERBLA.241"></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.244"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.245"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      DOUBLE PRECISION   <a name="DLAMCH.246"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANGE.246"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
      EXTERNAL           <a name="LSAME.247"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.247"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="DLAMCH.247"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANGE.247"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, INT, MAX
<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">     Decode the input arguments
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.256"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVL, <span class="string">'N'</span> ) ) THEN
         IJOBVL = 1
         ILVL = .FALSE.
      ELSE IF( <a name="LSAME.259"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVL, <span class="string">'V'</span> ) ) THEN
         IJOBVL = 2
         ILVL = .TRUE.
      ELSE
         IJOBVL = -1
         ILVL = .FALSE.
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.267"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVR, <span class="string">'N'</span> ) ) THEN
         IJOBVR = 1
         ILVR = .FALSE.
      ELSE IF( <a name="LSAME.270"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVR, <span class="string">'V'</span> ) ) THEN
         IJOBVR = 2
         ILVR = .TRUE.
      ELSE
         IJOBVR = -1
         ILVR = .FALSE.
      END IF
      ILV = ILVL .OR. ILVR
<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>      LWKMIN = MAX( 8*N, 1 )
      LWKOPT = LWKMIN
      WORK( 1 ) = LWKOPT
      LQUERY = ( LWORK.EQ.-1 )
      INFO = 0
      IF( IJOBVL.LE.0 ) THEN
         INFO = -1
      ELSE IF( IJOBVR.LE.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
         INFO = -12
      ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
         INFO = -14
      ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
         INFO = -16
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         NB1 = <a name="ILAENV.305"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DGEQRF.305"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>'</span>, <span class="string">' '</span>, N, N, -1, -1 )
         NB2 = <a name="ILAENV.306"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMQR.306"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>'</span>, <span class="string">' '</span>, N, N, N, -1 )
         NB3 = <a name="ILAENV.307"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORGQR.307"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>'</span>, <span class="string">' '</span>, N, N, N, -1 )
         NB = MAX( NB1, NB2, NB3 )
         LOPT = 2*N + MAX( 6*N, N*( NB+1 ) )
         WORK( 1 ) = LOPT
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.314"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DGEGV.314"></a><a href="dgegv.f.html#DGEGV.1">DGEGV</a> '</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.327"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'E'</span> )*<a name="DLAMCH.327"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'B'</span> )