zggevx.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) DOUBLE PRECISION 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) DOUBLE PRECISION
</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) DOUBLE PRECISION
</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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          If JOB = '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*16 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="ZHGEQZ.218"></a><a href="zhgeqz.f.html#ZHGEQZ.1">ZHGEQZ</a>.
</span><span class="comment">*</span><span class="comment">                =N+2: error return from <a name="ZTGEVC.219"></a><a href="ztgevc.f.html#ZTGEVC.1">ZTGEVC</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>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
      COMPLEX*16         CZERO, CONE
      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
     $                   CONE = ( 1.0D+0, 0.0D+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
      DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
     $                   SMLNUM, TEMP
      COMPLEX*16         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="DLABAD.268"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, <a name="DLASCL.268"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="XERBLA.268"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, <a name="ZGEQRF.268"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</a>, <a name="ZGGBAK.268"></a><a href="zggbak.f.html#ZGGBAK.1">ZGGBAK</a>, <a name="ZGGBAL.268"></a><a href="zggbal.f.html#ZGGBAL.1">ZGGBAL</a>,
     $                   <a name="ZGGHRD.269"></a><a href="zgghrd.f.html#ZGGHRD.1">ZGGHRD</a>, <a name="ZHGEQZ.269"></a><a href="zhgeqz.f.html#ZHGEQZ.1">ZHGEQZ</a>, <a name="ZLACPY.269"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLASCL.269"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>, <a name="ZLASET.269"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>, <a name="ZTGEVC.269"></a><a href="ztgevc.f.html#ZTGEVC.1">ZTGEVC</a>,
     $                   <a name="ZTGSNA.270"></a><a href="ztgsna.f.html#ZTGSNA.1">ZTGSNA</a>, <a name="ZUNGQR.270"></a><a href="zungqr.f.html#ZUNGQR.1">ZUNGQR</a>, <a name="ZUNMQR.270"></a><a href="zunmqr.f.html#ZUNMQR.1">ZUNMQR</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>
      DOUBLE PRECISION   <a name="DLAMCH.275"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGE.275"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</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="DLAMCH.276"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGE.276"></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          ABS, DBLE, DIMAG, MAX, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      DOUBLE PRECISION   ABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( 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">