slarre.f.html

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

<span class="comment">*</span><span class="comment">           check for element growth
</span>            IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
               NOREP = .TRUE.
            ELSE
               NOREP = .FALSE.
            ENDIF
            IF( USEDQD .AND. .NOT.NOREP ) THEN
<span class="comment">*</span><span class="comment">              Ensure the definiteness of the representation
</span><span class="comment">*</span><span class="comment">              All entries of D (of L D L^T) must have the same sign
</span>               DO 71 I = 1, IN
                  TMP = SGNDEF*WORK( I )
                  IF( TMP.LT.ZERO ) NOREP = .TRUE.
 71            CONTINUE
            ENDIF
            IF(NOREP) THEN
<span class="comment">*</span><span class="comment">              Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
</span><span class="comment">*</span><span class="comment">              shift which makes the matrix definite. So we should end up
</span><span class="comment">*</span><span class="comment">              here really only in the case of IRANGE = VALRNG or INDRNG.
</span>               IF( IDUM.EQ.MAXTRY-1 ) THEN
                  IF( SGNDEF.EQ.ONE ) THEN
<span class="comment">*</span><span class="comment">                    The fudged Gerschgorin shift should succeed
</span>                     SIGMA =
     $                    GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
                  ELSE
                     SIGMA =
     $                    GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
                  END IF
               ELSE
                  SIGMA = SIGMA - SGNDEF * TAU
                  TAU = TWO * TAU
               END IF
            ELSE
<span class="comment">*</span><span class="comment">              an initial RRR is found
</span>               GO TO 83
            END IF
 80      CONTINUE
<span class="comment">*</span><span class="comment">        if the program reaches this point, no base representation could be
</span><span class="comment">*</span><span class="comment">        found in MAXTRY iterations.
</span>         INFO = 2
         RETURN

 83      CONTINUE
<span class="comment">*</span><span class="comment">        At this point, we have found an initial base representation
</span><span class="comment">*</span><span class="comment">        T - SIGMA I = L D L^T with not too much element growth.
</span><span class="comment">*</span><span class="comment">        Store the shift.
</span>         E( IEND ) = SIGMA
<span class="comment">*</span><span class="comment">        Store D and L.
</span>         CALL SCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
         CALL SCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )


         IF(MB.GT.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Perturb each entry of the base representation by a small
</span><span class="comment">*</span><span class="comment">           (but random) relative amount to overcome difficulties with
</span><span class="comment">*</span><span class="comment">           glued matrices.
</span><span class="comment">*</span><span class="comment">
</span>            DO 122 I = 1, 4
               ISEED( I ) = 1
 122        CONTINUE

            CALL <a name="SLARNV.628"></a><a href="slarnv.f.html#SLARNV.1">SLARNV</a>(2, ISEED, 2*IN-1, WORK(1))
            DO 125 I = 1,IN-1
               D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
               E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
 125        CONTINUE
            D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
<span class="comment">*</span><span class="comment">
</span>         ENDIF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Don't update the Gerschgorin intervals because keeping track
</span><span class="comment">*</span><span class="comment">        of the updates would be too much work in <a name="SLARRV.638"></a><a href="slarrv.f.html#SLARRV.1">SLARRV</a>.
</span><span class="comment">*</span><span class="comment">        We update W instead and use it to locate the proper Gerschgorin
</span><span class="comment">*</span><span class="comment">        intervals.
</span>
<span class="comment">*</span><span class="comment">        Compute the required eigenvalues of L D L' by bisection or dqds
</span>         IF ( .NOT.USEDQD ) THEN
<span class="comment">*</span><span class="comment">           If <a name="SLARRD.644"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a> has been used, shift the eigenvalue approximations
</span><span class="comment">*</span><span class="comment">           according to their representation. This is necessary for
</span><span class="comment">*</span><span class="comment">           a uniform <a name="SLARRV.646"></a><a href="slarrv.f.html#SLARRV.1">SLARRV</a> since dqds computes eigenvalues of the
</span><span class="comment">*</span><span class="comment">           shifted representation. In <a name="SLARRV.647"></a><a href="slarrv.f.html#SLARRV.1">SLARRV</a>, W will always hold the
</span><span class="comment">*</span><span class="comment">           UNshifted eigenvalue approximation.
</span>            DO 134 J=WBEGIN,WEND
               W(J) = W(J) - SIGMA
               WERR(J) = WERR(J) + ABS(W(J)) * EPS
 134        CONTINUE
<span class="comment">*</span><span class="comment">           call <a name="SLARRB.653"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a> to reduce eigenvalue error of the approximations
</span><span class="comment">*</span><span class="comment">           from <a name="SLARRD.654"></a><a href="slarrd.f.html#SLARRD.1">SLARRD</a>
</span>            DO 135 I = IBEGIN, IEND-1
               WORK( I ) = D( I ) * E( I )**2
 135        CONTINUE
<span class="comment">*</span><span class="comment">           use bisection to find EV from INDL to INDU
</span>            CALL <a name="SLARRB.659"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a>(IN, D(IBEGIN), WORK(IBEGIN),
     $                  INDL, INDU, RTOL1, RTOL2, INDL-1,
     $                  W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
     $                  WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
     $                  IN, IINFO )
            IF( IINFO .NE. 0 ) THEN
               INFO = -4
               RETURN
            END IF
<span class="comment">*</span><span class="comment">           <a name="SLARRB.668"></a><a href="slarrb.f.html#SLARRB.1">SLARRB</a> computes all gaps correctly except for the last one
</span><span class="comment">*</span><span class="comment">           Record distance to VU/GU
</span>            WGAP( WEND ) = MAX( ZERO,
     $           ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
            DO 138 I = INDL, INDU
               M = M + 1
               IBLOCK(M) = JBLK
               INDEXW(M) = I
 138        CONTINUE
         ELSE
<span class="comment">*</span><span class="comment">           Call dqds to get all eigs (and then possibly delete unwanted
</span><span class="comment">*</span><span class="comment">           eigenvalues).
</span><span class="comment">*</span><span class="comment">           Note that dqds finds the eigenvalues of the L D L^T representation
</span><span class="comment">*</span><span class="comment">           of T to high relative accuracy. High relative accuracy
</span><span class="comment">*</span><span class="comment">           might be lost when the shift of the RRR is subtracted to obtain
</span><span class="comment">*</span><span class="comment">           the eigenvalues of T. However, T is not guaranteed to define its
</span><span class="comment">*</span><span class="comment">           eigenvalues to high relative accuracy anyway.
</span><span class="comment">*</span><span class="comment">           Set RTOL to the order of the tolerance used in <a name="SLASQ2.685"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>
</span><span class="comment">*</span><span class="comment">           This is an ESTIMATED error, the worst case bound is 4*N*EPS
</span><span class="comment">*</span><span class="comment">           which is usually too large and requires unnecessary work to be
</span><span class="comment">*</span><span class="comment">           done by bisection when computing the eigenvectors
</span>            RTOL = LOG(REAL(IN)) * FOUR * EPS
            J = IBEGIN
            DO 140 I = 1, IN - 1
               WORK( 2*I-1 ) = ABS( D( J ) )
               WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
               J = J + 1
  140       CONTINUE
            WORK( 2*IN-1 ) = ABS( D( IEND ) )
            WORK( 2*IN ) = ZERO
            CALL <a name="SLASQ2.698"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>( IN, WORK, IINFO )
            IF( IINFO .NE. 0 ) THEN
<span class="comment">*</span><span class="comment">              If IINFO = -5 then an index is part of a tight cluster
</span><span class="comment">*</span><span class="comment">              and should be changed. The index is in IWORK(1) and the
</span><span class="comment">*</span><span class="comment">              gap is in WORK(N+1)
</span>               INFO = -5
               RETURN
            ELSE
<span class="comment">*</span><span class="comment">              Test that all eigenvalues are positive as expected
</span>               DO 149 I = 1, IN
                  IF( WORK( I ).LT.ZERO ) THEN
                     INFO = -6
                     RETURN
                  ENDIF
 149           CONTINUE
            END IF
            IF( SGNDEF.GT.ZERO ) THEN
               DO 150 I = INDL, INDU
                  M = M + 1
                  W( M ) = WORK( IN-I+1 )
                  IBLOCK( M ) = JBLK
                  INDEXW( M ) = I
 150           CONTINUE
            ELSE
               DO 160 I = INDL, INDU
                  M = M + 1
                  W( M ) = -WORK( I )
                  IBLOCK( M ) = JBLK
                  INDEXW( M ) = I
 160           CONTINUE
            END IF

            DO 165 I = M - MB + 1, M
<span class="comment">*</span><span class="comment">              the value of RTOL below should be the tolerance in <a name="SLASQ2.731"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>
</span>               WERR( I ) = RTOL * ABS( W(I) )
 165        CONTINUE
            DO 166 I = M - MB + 1, M - 1
<span class="comment">*</span><span class="comment">              compute the right gap between the intervals
</span>               WGAP( I ) = MAX( ZERO,
     $                          W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
 166        CONTINUE
            WGAP( M ) = MAX( ZERO,
     $           ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
         END IF
<span class="comment">*</span><span class="comment">        proceed with next block
</span>         IBEGIN = IEND + 1
         WBEGIN = WEND + 1
 170  CONTINUE
<span class="comment">*</span><span class="comment">
</span>
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     end of <a name="SLARRE.750"></a><a href="slarre.f.html#SLARRE.1">SLARRE</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>