sgegv.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="SGEQRF.160"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>, <a name="SORMQR.160"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>, and <a name="SORGQR.160"></a><a href="sorgqr.f.html#SORGQR.1">SORGQR</a>.)  Then compute:
</span><span class="comment">*</span><span class="comment">          NB  -- MAX of the blocksizes for <a name="SGEQRF.161"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>, <a name="SORMQR.161"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>, and <a name="SORGQR.161"></a><a href="sorgqr.f.html#SORGQR.1">SORGQR</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="SGGBAL.178"></a><a href="sggbal.f.html#SGGBAL.1">SGGBAL</a>
</span><span class="comment">*</span><span class="comment">                =N+2: error return from <a name="SGEQRF.179"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>
</span><span class="comment">*</span><span class="comment">                =N+3: error return from <a name="SORMQR.180"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>
</span><span class="comment">*</span><span class="comment">                =N+4: error return from <a name="SORGQR.181"></a><a href="sorgqr.f.html#SORGQR.1">SORGQR</a>
</span><span class="comment">*</span><span class="comment">                =N+5: error return from <a name="SGGHRD.182"></a><a href="sgghrd.f.html#SGGHRD.1">SGGHRD</a>
</span><span class="comment">*</span><span class="comment">                =N+6: error return from <a name="SHGEQZ.183"></a><a href="shgeqz.f.html#SHGEQZ.1">SHGEQZ</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="STGEVC.185"></a><a href="stgevc.f.html#STGEVC.1">STGEVC</a>
</span><span class="comment">*</span><span class="comment">                =N+8: error return from <a name="SGGBAK.186"></a><a href="sggbak.f.html#SGGBAK.1">SGGBAK</a> (computing VL)
</span><span class="comment">*</span><span class="comment">                =N+9: error return from <a name="SGGBAK.187"></a><a href="sggbak.f.html#SGGBAK.1">SGGBAK</a> (computing VR)
</span><span class="comment">*</span><span class="comment">                =N+10: error return from <a name="SLASCL.188"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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="SGGBAL.196"></a><a href="sggbal.f.html#SGGBAL.1">SGGBAL</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="SGGBAK.205"></a><a href="sggbak.f.html#SGGBAK.1">SGGBAK</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="SHGEQZ.217"></a><a href="shgeqz.f.html#SHGEQZ.1">SHGEQZ</a>, <a name="SGEGS.217"></a><a href="sgegs.f.html#SGEGS.1">SGEGS</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>      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            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
      REAL               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="SGEQRF.240"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>, <a name="SGGBAK.240"></a><a href="sggbak.f.html#SGGBAK.1">SGGBAK</a>, <a name="SGGBAL.240"></a><a href="sggbal.f.html#SGGBAL.1">SGGBAL</a>, <a name="SGGHRD.240"></a><a href="sgghrd.f.html#SGGHRD.1">SGGHRD</a>, <a name="SHGEQZ.240"></a><a href="shgeqz.f.html#SHGEQZ.1">SHGEQZ</a>, <a name="SLACPY.240"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>,
     $                   <a name="SLASCL.241"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>, <a name="SLASET.241"></a><a href="slaset.f.html#SLASET.1">SLASET</a>, <a name="SORGQR.241"></a><a href="sorgqr.f.html#SORGQR.1">SORGQR</a>, <a name="SORMQR.241"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>, <a name="STGEVC.241"></a><a href="stgevc.f.html#STGEVC.1">STGEVC</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>
      REAL               <a name="SLAMCH.246"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANGE.246"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>
      EXTERNAL           <a name="ILAENV.247"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="LSAME.247"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="SLAMCH.247"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANGE.247"></a><a href="slange.f.html#SLANGE.1">SLANGE</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="SGEQRF.305"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</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="SORMQR.306"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</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="SORGQR.307"></a><a href="sorgqr.f.html#SORGQR.1">SORGQR</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="SGEGV.314"></a><a href="sgegv.f.html#SGEGV.1">SGEGV</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="SLAMCH.327"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'E'</span> )*<a name="SLAMCH.327"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'B'</span> )