dhgeqz.f.html

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

     $   GO TO 380
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     MAIN QZ ITERATION LOOP
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize dynamic indices
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Eigenvalues ILAST+1:N have been found.
</span><span class="comment">*</span><span class="comment">        Column operations modify rows IFRSTM:whatever.
</span><span class="comment">*</span><span class="comment">        Row operations modify columns whatever:ILASTM.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If only eigenvalues are being computed, then
</span><span class="comment">*</span><span class="comment">        IFRSTM is the row of the last splitting row above row ILAST;
</span><span class="comment">*</span><span class="comment">        this is always at least ILO.
</span><span class="comment">*</span><span class="comment">     IITER counts iterations since the last eigenvalue was found,
</span><span class="comment">*</span><span class="comment">        to tell when to use an extraordinary shift.
</span><span class="comment">*</span><span class="comment">     MAXIT is the maximum number of QZ sweeps allowed.
</span><span class="comment">*</span><span class="comment">
</span>      ILAST = IHI
      IF( ILSCHR ) THEN
         IFRSTM = 1
         ILASTM = N
      ELSE
         IFRSTM = ILO
         ILASTM = IHI
      END IF
      IITER = 0
      ESHIFT = ZERO
      MAXIT = 30*( IHI-ILO+1 )
<span class="comment">*</span><span class="comment">
</span>      DO 360 JITER = 1, MAXIT
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Split the matrix if possible.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Two tests:
</span><span class="comment">*</span><span class="comment">           1: H(j,j-1)=0  or  j=ILO
</span><span class="comment">*</span><span class="comment">           2: T(j,j)=0
</span><span class="comment">*</span><span class="comment">
</span>         IF( ILAST.EQ.ILO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Special case: j=ILAST
</span><span class="comment">*</span><span class="comment">
</span>            GO TO 80
         ELSE
            IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
               H( ILAST, ILAST-1 ) = ZERO
               GO TO 80
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
            T( ILAST, ILAST ) = ZERO
            GO TO 70
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        General case: j&lt;ILAST
</span><span class="comment">*</span><span class="comment">
</span>         DO 60 J = ILAST - 1, ILO, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test 1: for H(j,j-1)=0 or j=ILO
</span><span class="comment">*</span><span class="comment">
</span>            IF( J.EQ.ILO ) THEN
               ILAZRO = .TRUE.
            ELSE
               IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN
                  H( J, J-1 ) = ZERO
                  ILAZRO = .TRUE.
               ELSE
                  ILAZRO = .FALSE.
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test 2: for T(j,j)=0
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( T( J, J ) ).LT.BTOL ) THEN
               T( J, J ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Test 1a: Check for 2 consecutive small subdiagonals in A
</span><span class="comment">*</span><span class="comment">
</span>               ILAZR2 = .FALSE.
               IF( .NOT.ILAZRO ) THEN
                  TEMP = ABS( H( J, J-1 ) )
                  TEMP2 = ABS( H( J, J ) )
                  TEMPR = MAX( TEMP, TEMP2 )
                  IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
                     TEMP = TEMP / TEMPR
                     TEMP2 = TEMP2 / TEMPR
                  END IF
                  IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2*
     $                ( ASCALE*ATOL ) )ILAZR2 = .TRUE.
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              If both tests pass (1 &amp; 2), i.e., the leading diagonal
</span><span class="comment">*</span><span class="comment">              element of B in the block is zero, split a 1x1 block off
</span><span class="comment">*</span><span class="comment">              at the top. (I.e., at the J-th row/column) The leading
</span><span class="comment">*</span><span class="comment">              diagonal element of the remainder can also be zero, so
</span><span class="comment">*</span><span class="comment">              this may have to be done repeatedly.
</span><span class="comment">*</span><span class="comment">
</span>               IF( ILAZRO .OR. ILAZR2 ) THEN
                  DO 40 JCH = J, ILAST - 1
                     TEMP = H( JCH, JCH )
                     CALL <a name="DLARTG.478"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( TEMP, H( JCH+1, JCH ), C, S,
     $                            H( JCH, JCH ) )
                     H( JCH+1, JCH ) = ZERO
                     CALL DROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
     $                          H( JCH+1, JCH+1 ), LDH, C, S )
                     CALL DROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
     $                          T( JCH+1, JCH+1 ), LDT, C, S )
                     IF( ILQ )
     $                  CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
     $                             C, S )
                     IF( ILAZR2 )
     $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
                     ILAZR2 = .FALSE.
                     IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
                        IF( JCH+1.GE.ILAST ) THEN
                           GO TO 80
                        ELSE
                           IFIRST = JCH + 1
                           GO TO 110
                        END IF
                     END IF
                     T( JCH+1, JCH+1 ) = ZERO
   40             CONTINUE
                  GO TO 70
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
</span><span class="comment">*</span><span class="comment">                 Then process as in the case T(ILAST,ILAST)=0
</span><span class="comment">*</span><span class="comment">
</span>                  DO 50 JCH = J, ILAST - 1
                     TEMP = T( JCH, JCH+1 )
                     CALL <a name="DLARTG.509"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( TEMP, T( JCH+1, JCH+1 ), C, S,
     $                            T( JCH, JCH+1 ) )
                     T( JCH+1, JCH+1 ) = ZERO
                     IF( JCH.LT.ILASTM-1 )
     $                  CALL DROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
     $                             T( JCH+1, JCH+2 ), LDT, C, S )
                     CALL DROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
     $                          H( JCH+1, JCH-1 ), LDH, C, S )
                     IF( ILQ )
     $                  CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
     $                             C, S )
                     TEMP = H( JCH+1, JCH )
                     CALL <a name="DLARTG.521"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( TEMP, H( JCH+1, JCH-1 ), C, S,
     $                            H( JCH+1, JCH ) )
                     H( JCH+1, JCH-1 ) = ZERO
                     CALL DROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
     $                          H( IFRSTM, JCH-1 ), 1, C, S )
                     CALL DROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
     $                          T( IFRSTM, JCH-1 ), 1, C, S )
                     IF( ILZ )
     $                  CALL DROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
     $                             C, S )
   50             CONTINUE
                  GO TO 70
               END IF
            ELSE IF( ILAZRO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Only test 1 passed -- work on J:ILAST
</span><span class="comment">*</span><span class="comment">
</span>               IFIRST = J
               GO TO 110
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Neither test passed -- try next J
</span><span class="comment">*</span><span class="comment">
</span>   60    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        (Drop-through is &quot;impossible&quot;)
</span><span class="comment">*</span><span class="comment">
</span>         INFO = N + 1
         GO TO 420
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
</span><span class="comment">*</span><span class="comment">        1x1 block.
</span><span class="comment">*</span><span class="comment">
</span>   70    CONTINUE
         TEMP = H( ILAST, ILAST )
         CALL <a name="DLARTG.556"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( TEMP, H( ILAST, ILAST-1 ), C, S,
     $                H( ILAST, ILAST ) )
         H( ILAST, ILAST-1 ) = ZERO
         CALL DROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
     $              H( IFRSTM, ILAST-1 ), 1, C, S )
         CALL DROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
     $              T( IFRSTM, ILAST-1 ), 1, C, S )
         IF( ILZ )
     $      CALL DROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
</span><span class="comment">*</span><span class="comment">                              and BETA
</span><span class="comment">*</span><span class="comment">
</span>   80    CONTINUE
         IF( T( ILAST, ILAST ).LT.ZERO ) THEN
            IF( ILSCHR ) THEN
               DO 90 J = IFRSTM, ILAST
                  H( J, ILAST ) = -H( J, ILAST )
                  T( J, ILAST ) = -T( J, ILAST )
   90          CONTINUE