cggevx.f.html

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

</span><span class="comment">*</span><span class="comment">                      = DL(j)  for j = ILO,...,IHI
</span><span class="comment">*</span><span class="comment">                      = PL(j)  for j = IHI+1,...,N.
</span><span class="comment">*</span><span class="comment">          The order in which the interchanges are made is N to IHI+1,
</span><span class="comment">*</span><span class="comment">          then 1 to ILO-1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RSCALE  (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          Details of the permutations and scaling factors applied
</span><span class="comment">*</span><span class="comment">          to the right side of A and B.  If PR(j) is the index of the
</span><span class="comment">*</span><span class="comment">          column interchanged with column j, and DR(j) is the scaling
</span><span class="comment">*</span><span class="comment">          factor applied to column j, then
</span><span class="comment">*</span><span class="comment">            RSCALE(j) = PR(j)  for j = 1,...,ILO-1
</span><span class="comment">*</span><span class="comment">                      = DR(j)  for j = ILO,...,IHI
</span><span class="comment">*</span><span class="comment">                      = PR(j)  for j = IHI+1,...,N
</span><span class="comment">*</span><span class="comment">          The order in which the interchanges are made is N to IHI+1,
</span><span class="comment">*</span><span class="comment">          then 1 to ILO-1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ABNRM   (output) REAL
</span><span class="comment">*</span><span class="comment">          The one-norm of the balanced matrix A.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  BBNRM   (output) REAL
</span><span class="comment">*</span><span class="comment">          The one-norm of the balanced matrix B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RCONDE  (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          If SENSE = 'E' or 'B', the reciprocal condition numbers of
</span><span class="comment">*</span><span class="comment">          the eigenvalues, stored in consecutive elements of the array.
</span><span class="comment">*</span><span class="comment">          If SENSE = 'N' or 'V', RCONDE is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RCONDV  (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          If SENSE = 'V' or 'B', the estimated reciprocal condition
</span><span class="comment">*</span><span class="comment">          numbers of the eigenvectors, stored in consecutive elements
</span><span class="comment">*</span><span class="comment">          of the array. If the eigenvalues cannot be reordered to
</span><span class="comment">*</span><span class="comment">          compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
</span><span class="comment">*</span><span class="comment">          when the true value would be very small anyway. 
</span><span class="comment">*</span><span class="comment">          If SENSE = 'N' or 'E', RCONDV is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
</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,2*N).
</span><span class="comment">*</span><span class="comment">          If SENSE = 'E', LWORK &gt;= max(1,4*N).
</span><span class="comment">*</span><span class="comment">          If SENSE = 'V' or 'B', LWORK &gt;= max(1,2*N*N+2*N).
</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.198"></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 (lrwork)
</span><span class="comment">*</span><span class="comment">          lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
</span><span class="comment">*</span><span class="comment">          and at least max(1,2*N) otherwise.
</span><span class="comment">*</span><span class="comment">          Real workspace.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace) INTEGER array, dimension (N+2)
</span><span class="comment">*</span><span class="comment">          If SENSE = 'E', IWORK is not referenced.
</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">          If SENSE = 'N', BWORK is not referenced.
</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 ALPHA(j) and BETA(j) should be correct
</span><span class="comment">*</span><span class="comment">                for j=INFO+1,...,N.
</span><span class="comment">*</span><span class="comment">          &gt; N:  =N+1: other than QZ iteration failed in <a name="CHGEQZ.218"></a><a href="chgeqz.f.html#CHGEQZ.1">CHGEQZ</a>.
</span><span class="comment">*</span><span class="comment">                =N+2: error return from <a name="CTGEVC.219"></a><a href="ctgevc.f.html#CTGEVC.1">CTGEVC</a>.
</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 a matrix pair (A,B) includes, first, permuting rows and
</span><span class="comment">*</span><span class="comment">  columns to isolate eigenvalues, second, applying diagonal similarity
</span><span class="comment">*</span><span class="comment">  transformation to the rows and columns to make the rows and columns
</span><span class="comment">*</span><span class="comment">  as close in norm as possible. The computed reciprocal condition
</span><span class="comment">*</span><span class="comment">  numbers correspond to the balanced matrix. Permuting rows and columns
</span><span class="comment">*</span><span class="comment">  will not change the condition numbers (in exact arithmetic) but
</span><span class="comment">*</span><span class="comment">  diagonal scaling will.  For further explanation of balancing, see
</span><span class="comment">*</span><span class="comment">  section 4.11.1.2 of LAPACK Users' Guide.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate error bound on the chordal distance between the i-th
</span><span class="comment">*</span><span class="comment">  computed generalized eigenvalue w and the corresponding exact
</span><span class="comment">*</span><span class="comment">  eigenvalue lambda is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       chord(w, lambda) &lt;= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate error bound for the angle between the i-th computed
</span><span class="comment">*</span><span class="comment">  eigenvector VL(i) or VR(i) is given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       EPS * norm(ABNRM, BBNRM) / DIF(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  For further explanation of the reciprocal condition numbers RCONDE
</span><span class="comment">*</span><span class="comment">  and RCONDV, see section 4.11 of LAPACK User's Guide.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
      COMPLEX            CZERO, CONE
      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
     $                   CONE = ( 1.0E+0, 0.0E+0 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
     $                   WANTSB, WANTSE, WANTSN, WANTSV
      CHARACTER          CHTEMP
      INTEGER            I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
     $                   ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK
      REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
     $                   SMLNUM, TEMP
      COMPLEX            X
<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="CGEQRF.268"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</a>, <a name="CGGBAK.268"></a><a href="cggbak.f.html#CGGBAK.1">CGGBAK</a>, <a name="CGGBAL.268"></a><a href="cggbal.f.html#CGGBAL.1">CGGBAL</a>, <a name="CGGHRD.268"></a><a href="cgghrd.f.html#CGGHRD.1">CGGHRD</a>, <a name="CHGEQZ.268"></a><a href="chgeqz.f.html#CHGEQZ.1">CHGEQZ</a>, <a name="CLACPY.268"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>,
     $                   <a name="CLASCL.269"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>, <a name="CLASET.269"></a><a href="claset.f.html#CLASET.1">CLASET</a>, <a name="CTGEVC.269"></a><a href="ctgevc.f.html#CTGEVC.1">CTGEVC</a>, <a name="CTGSNA.269"></a><a href="ctgsna.f.html#CTGSNA.1">CTGSNA</a>, <a name="CUNGQR.269"></a><a href="cungqr.f.html#CUNGQR.1">CUNGQR</a>, <a name="CUNMQR.269"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>,
     $                   <a name="SLABAD.270"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>, <a name="SLASCL.270"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>, <a name="XERBLA.270"></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.273"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.274"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      REAL               <a name="CLANGE.275"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>, <a name="SLAMCH.275"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="LSAME.276"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.276"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="CLANGE.276"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>, <a name="SLAMCH.276"></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          ABS, AIMAG, MAX, REAL, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      REAL               ABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
<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.291"></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.294"></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.302"></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.305"></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>      NOSCL  = <a name="LSAME.314"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( BALANC, <span class="string">'N'</span> ) .OR. <a name="LSAME.314"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( BALANC, <span class="string">'P'</span> )
      WANTSN = <a name="LSAME.315"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SENSE, <span class="string">'N'</span> )
      WANTSE = <a name="LSAME.316"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SENSE, <span class="string">'E'</span> )
      WANTSV = <a name="LSAME.317"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SENSE, <span class="string">'V'</span> )
      WANTSB = <a name="LSAME.318"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SENSE, <span class="string">'B'</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">