claqr3.f.html

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

     $        THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== One more converged eigenvalue ====
</span><span class="comment">*</span><span class="comment">
</span>            NS = NS - 1
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== One undflatable eigenvalue.  Move it up out of the
</span><span class="comment">*</span><span class="comment">           .    way.   (<a name="CTREXC.306"></a><a href="ctrexc.f.html#CTREXC.1">CTREXC</a> can not fail in this case.) ====
</span><span class="comment">*</span><span class="comment">
</span>            IFST = NS
            CALL <a name="CTREXC.309"></a><a href="ctrexc.f.html#CTREXC.1">CTREXC</a>( <span class="string">'V'</span>, JW, T, LDT, V, LDV, IFST, ILST, INFO )
            ILST = ILST + 1
         END IF
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Return to Hessenberg form ====
</span><span class="comment">*</span><span class="comment">
</span>      IF( NS.EQ.0 )
     $   S = ZERO
<span class="comment">*</span><span class="comment">
</span>      IF( NS.LT.JW ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== sorting the diagonal of T improves accuracy for
</span><span class="comment">*</span><span class="comment">        .    graded matrices.  ====
</span><span class="comment">*</span><span class="comment">
</span>         DO 30 I = INFQR + 1, NS
            IFST = I
            DO 20 J = I + 1, NS
               IF( CABS1( T( J, J ) ).GT.CABS1( T( IFST, IFST ) ) )
     $            IFST = J
   20       CONTINUE
            ILST = I
            IF( IFST.NE.ILST )
     $         CALL <a name="CTREXC.332"></a><a href="ctrexc.f.html#CTREXC.1">CTREXC</a>( <span class="string">'V'</span>, JW, T, LDT, V, LDV, IFST, ILST, INFO )
   30    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Restore shift/eigenvalue array from T ====
</span><span class="comment">*</span><span class="comment">
</span>      DO 40 I = INFQR + 1, JW
         SH( KWTOP+I-1 ) = T( I, I )
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>      IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
         IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Reflect spike back into lower triangle ====
</span><span class="comment">*</span><span class="comment">
</span>            CALL CCOPY( NS, V, LDV, WORK, 1 )
            DO 50 I = 1, NS
               WORK( I ) = CONJG( WORK( I ) )
   50       CONTINUE
            BETA = WORK( 1 )
            CALL <a name="CLARFG.353"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( NS, BETA, WORK( 2 ), 1, TAU )
            WORK( 1 ) = ONE
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="CLASET.356"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'L'</span>, JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="CLARF.358"></a><a href="clarf.f.html#CLARF.1">CLARF</a>( <span class="string">'L'</span>, NS, JW, WORK, 1, CONJG( TAU ), T, LDT,
     $                  WORK( JW+1 ) )
            CALL <a name="CLARF.360"></a><a href="clarf.f.html#CLARF.1">CLARF</a>( <span class="string">'R'</span>, NS, NS, WORK, 1, TAU, T, LDT,
     $                  WORK( JW+1 ) )
            CALL <a name="CLARF.362"></a><a href="clarf.f.html#CLARF.1">CLARF</a>( <span class="string">'R'</span>, JW, NS, WORK, 1, TAU, V, LDV,
     $                  WORK( JW+1 ) )
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="CGEHRD.365"></a><a href="cgehrd.f.html#CGEHRD.1">CGEHRD</a>( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
     $                   LWORK-JW, INFO )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Copy updated reduced window into place ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( KWTOP.GT.1 )
     $      H( KWTOP, KWTOP-1 ) = S*CONJG( V( 1, 1 ) )
         CALL <a name="CLACPY.373"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'U'</span>, JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )
         CALL CCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
     $               LDH+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Accumulate orthogonal matrix in order update
</span><span class="comment">*</span><span class="comment">        .    H and Z, if requested.  (A modified version
</span><span class="comment">*</span><span class="comment">        .    of  <a name="CUNGHR.379"></a><a href="cunghr.f.html#CUNGHR.1">CUNGHR</a> that accumulates block Householder
</span><span class="comment">*</span><span class="comment">        .    transformations into V directly might be
</span><span class="comment">*</span><span class="comment">        .    marginally more efficient than the following.) ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
            CALL <a name="CUNGHR.384"></a><a href="cunghr.f.html#CUNGHR.1">CUNGHR</a>( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
     $                   LWORK-JW, INFO )
            CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JW, NS, NS, ONE, V, LDV, T, LDT, ZERO,
     $                  WV, LDWV )
            CALL <a name="CLACPY.388"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'A'</span>, JW, NS, WV, LDWV, V, LDV )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Update vertical slab in H ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( WANTT ) THEN
            LTOP = 1
         ELSE
            LTOP = KTOP
         END IF
         DO 60 KROW = LTOP, KWTOP - 1, NV
            KLN = MIN( NV, KWTOP-KROW )
            CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, KLN, JW, JW, ONE, H( KROW, KWTOP ),
     $                  LDH, V, LDV, ZERO, WV, LDWV )
            CALL <a name="CLACPY.402"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'A'</span>, KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )
   60    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Update horizontal slab in H ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( WANTT ) THEN
            DO 70 KCOL = KBOT + 1, N, NH
               KLN = MIN( NH, N-KCOL+1 )
               CALL CGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, JW, KLN, JW, ONE, V, LDV,
     $                     H( KWTOP, KCOL ), LDH, ZERO, T, LDT )
               CALL <a name="CLACPY.412"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'A'</span>, JW, KLN, T, LDT, H( KWTOP, KCOL ),
     $                      LDH )
   70       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Update vertical slab in Z ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( WANTZ ) THEN
            DO 80 KROW = ILOZ, IHIZ, NV
               KLN = MIN( NV, IHIZ-KROW+1 )
               CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, KLN, JW, JW, ONE, Z( KROW, KWTOP ),
     $                     LDZ, V, LDV, ZERO, WV, LDWV )
               CALL <a name="CLACPY.424"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'A'</span>, KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
     $                      LDZ )
   80       CONTINUE
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Return the number of deflations ... ====
</span><span class="comment">*</span><span class="comment">
</span>      ND = JW - NS
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== ... and the number of shifts. (Subtracting
</span><span class="comment">*</span><span class="comment">     .    INFQR from the spike length takes care
</span><span class="comment">*</span><span class="comment">     .    of the case of a rare QR failure while
</span><span class="comment">*</span><span class="comment">     .    calculating eigenvalues of the deflation
</span><span class="comment">*</span><span class="comment">     .    window.)  ====
</span><span class="comment">*</span><span class="comment">
</span>      NS = NS - INFQR
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      ==== Return optimal workspace. ====
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = CMPLX( LWKOPT, 0 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== End of <a name="CLAQR3.446"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</a> ====
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>