dtgsen.f.html

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

</span>      DOUBLE PRECISION   <a name="DLAMCH.356"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DLAMCH.357"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX, SIGN, 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">     Decode and test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
<span class="comment">*</span><span class="comment">
</span>      IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -9
      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
         INFO = -14
      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
         INFO = -16
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.384"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DTGSEN.384"></a><a href="dtgsen.f.html#DTGSEN.1">DTGSEN</a>'</span>, -INFO )
         RETURN
      END IF
<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.390"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="DLAMCH.391"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> ) / EPS
      IERR = 0
<span class="comment">*</span><span class="comment">
</span>      WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
      WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
      WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
      WANTD = WANTD1 .OR. WANTD2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set M to the dimension of the specified pair of deflating
</span><span class="comment">*</span><span class="comment">     subspaces.
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      PAIR = .FALSE.
      DO 10 K = 1, N
         IF( PAIR ) THEN
            PAIR = .FALSE.
         ELSE
            IF( K.LT.N ) THEN
               IF( A( K+1, K ).EQ.ZERO ) THEN
                  IF( SELECT( K ) )
     $               M = M + 1
               ELSE
                  PAIR = .TRUE.
                  IF( SELECT( K ) .OR. SELECT( K+1 ) )
     $               M = M + 2
               END IF
            ELSE
               IF( SELECT( N ) )
     $            M = M + 1
            END IF
         END IF
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
         LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) )
         LIWMIN = MAX( 1, N+6 )
      ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
         LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) )
         LIWMIN = MAX( 1, 2*M*( N-M ), N+6 )
      ELSE
         LWMIN = MAX( 1, 4*N+16 )
         LIWMIN = 1
      END IF
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWMIN
      IWORK( 1 ) = LIWMIN
<span class="comment">*</span><span class="comment">
</span>      IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -22
      ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -24
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.445"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DTGSEN.445"></a><a href="dtgsen.f.html#DTGSEN.1">DTGSEN</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( M.EQ.N .OR. M.EQ.0 ) THEN
         IF( WANTP ) THEN
            PL = ONE
            PR = ONE
         END IF
         IF( WANTD ) THEN
            DSCALE = ZERO
            DSUM = ONE
            DO 20 I = 1, N
               CALL <a name="DLASSQ.462"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>( N, A( 1, I ), 1, DSCALE, DSUM )
               CALL <a name="DLASSQ.463"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>( N, B( 1, I ), 1, DSCALE, DSUM )
   20       CONTINUE
            DIF( 1 ) = DSCALE*SQRT( DSUM )
            DIF( 2 ) = DIF( 1 )
         END IF
         GO TO 60
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Collect the selected blocks at the top-left corner of (A, B).
</span><span class="comment">*</span><span class="comment">
</span>      KS = 0
      PAIR = .FALSE.
      DO 30 K = 1, N
         IF( PAIR ) THEN
            PAIR = .FALSE.
         ELSE
<span class="comment">*</span><span class="comment">
</span>            SWAP = SELECT( K )
            IF( K.LT.N ) THEN
               IF( A( K+1, K ).NE.ZERO ) THEN
                  PAIR = .TRUE.
                  SWAP = SWAP .OR. SELECT( K+1 )
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( SWAP ) THEN
               KS = KS + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Swap the K-th block to position KS.
</span><span class="comment">*</span><span class="comment">              Perform the reordering of diagonal blocks in (A, B)
</span><span class="comment">*</span><span class="comment">              by orthogonal transformation matrices and update
</span><span class="comment">*</span><span class="comment">              Q and Z accordingly (if requested):
</span><span class="comment">*</span><span class="comment">
</span>               KK = K
               IF( K.NE.KS )
     $            CALL <a name="DTGEXC.498"></a><a href="dtgexc.f.html#DTGEXC.1">DTGEXC</a>( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
     $                         Z, LDZ, KK, KS, WORK, LWORK, IERR )
<span class="comment">*</span><span class="comment">
</span>               IF( IERR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Swap is rejected: exit.
</span><span class="comment">*</span><span class="comment">
</span>                  INFO = 1
                  IF( WANTP ) THEN
                     PL = ZERO
                     PR = ZERO
                  END IF
                  IF( WANTD ) THEN
                     DIF( 1 ) = ZERO
                     DIF( 2 ) = ZERO
                  END IF
                  GO TO 60
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( PAIR )
     $            KS = KS + 1
            END IF
         END IF
   30 CONTINUE
      IF( WANTP ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve generalized Sylvester equation for R and L
</span><span class="comment">*</span><span class="comment">        and compute PL and PR.
</span><span class="comment">*</span><span class="comment">
</span>         N1 = M
         N2 = N - M
         I = N1 + 1
         IJB = 0
         CALL <a name="DLACPY.531"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, N1, N2, A( 1, I ), LDA, WORK, N1 )
         CALL <a name="DLACPY.532"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
     $                N1 )
         CALL <a name="DTGSYL.534"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</a>( <span class="string">'N'</span>, IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
     $                N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
     $                DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
     $                LWORK-2*N1*N2, IWORK, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Estimate the reciprocal of norms of &quot;projections&quot; onto left
</span><span class="comment">*</span><span class="comment">        and right eigenspaces.
</span><span class="comment">*</span><span class="comment">
</span>         RDSCAL = ZERO