dlarrf.f.html

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

      SAWNAN1 = .FALSE.
      SAWNAN2 = .FALSE.
<span class="comment">*</span><span class="comment">     Ensure that we do not back off too much of the initial shifts
</span>      LDELTA = MIN(LDMAX,LDELTA)
      RDELTA = MIN(RDMAX,RDELTA)

<span class="comment">*</span><span class="comment">     Compute the element growth when shifting to both ends of the cluster
</span><span class="comment">*</span><span class="comment">     accept the shift if there is no element growth at one of the two ends
</span>
<span class="comment">*</span><span class="comment">     Left end
</span>      S = -LSIGMA
      DPLUS( 1 ) = D( 1 ) + S
      IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
         DPLUS(1) = -PIVMIN
<span class="comment">*</span><span class="comment">        Need to set SAWNAN1 because refined RRR test should not be used
</span><span class="comment">*</span><span class="comment">        in this case
</span>         SAWNAN1 = .TRUE.
      ENDIF
      MAX1 = ABS( DPLUS( 1 ) )
      DO 6 I = 1, N - 1
         LPLUS( I ) = LD( I ) / DPLUS( I )
         S = S*LPLUS( I )*L( I ) - LSIGMA
         DPLUS( I+1 ) = D( I+1 ) + S
         IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
            DPLUS(I+1) = -PIVMIN
<span class="comment">*</span><span class="comment">           Need to set SAWNAN1 because refined RRR test should not be used
</span><span class="comment">*</span><span class="comment">           in this case
</span>            SAWNAN1 = .TRUE.
         ENDIF
         MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
 6    CONTINUE
      SAWNAN1 = SAWNAN1 .OR.  <a name="DISNAN.212"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>( MAX1 )

      IF( FORCER .OR.
     $   (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
         SIGMA = LSIGMA
         SHIFT = SLEFT
         GOTO 100
      ENDIF

<span class="comment">*</span><span class="comment">     Right end
</span>      S = -RSIGMA
      WORK( 1 ) = D( 1 ) + S
      IF(ABS(WORK(1)).LT.PIVMIN) THEN
         WORK(1) = -PIVMIN
<span class="comment">*</span><span class="comment">        Need to set SAWNAN2 because refined RRR test should not be used
</span><span class="comment">*</span><span class="comment">        in this case
</span>         SAWNAN2 = .TRUE.
      ENDIF
      MAX2 = ABS( WORK( 1 ) )
      DO 7 I = 1, N - 1
         WORK( N+I ) = LD( I ) / WORK( I )
         S = S*WORK( N+I )*L( I ) - RSIGMA
         WORK( I+1 ) = D( I+1 ) + S
         IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
            WORK(I+1) = -PIVMIN
<span class="comment">*</span><span class="comment">           Need to set SAWNAN2 because refined RRR test should not be used
</span><span class="comment">*</span><span class="comment">           in this case
</span>            SAWNAN2 = .TRUE.
         ENDIF
         MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
 7    CONTINUE
      SAWNAN2 = SAWNAN2 .OR.  <a name="DISNAN.243"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>( MAX2 )

      IF( FORCER .OR.
     $   (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
         SIGMA = RSIGMA
         SHIFT = SRIGHT
         GOTO 100
      ENDIF
<span class="comment">*</span><span class="comment">     If we are at this point, both shifts led to too much element growth
</span>
<span class="comment">*</span><span class="comment">     Record the better of the two shifts (provided it didn't lead to NaN)
</span>      IF(SAWNAN1.AND.SAWNAN2) THEN
<span class="comment">*</span><span class="comment">        both MAX1 and MAX2 are NaN
</span>         GOTO 50
      ELSE
         IF( .NOT.SAWNAN1 ) THEN
            INDX = 1
            IF(MAX1.LE.SMLGROWTH) THEN
               SMLGROWTH = MAX1
               BESTSHIFT = LSIGMA
            ENDIF
         ENDIF
         IF( .NOT.SAWNAN2 ) THEN
            IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
            IF(MAX2.LE.SMLGROWTH) THEN
               SMLGROWTH = MAX2
               BESTSHIFT = RSIGMA
            ENDIF
         ENDIF
      ENDIF

<span class="comment">*</span><span class="comment">     If we are here, both the left and the right shift led to
</span><span class="comment">*</span><span class="comment">     element growth. If the element growth is moderate, then
</span><span class="comment">*</span><span class="comment">     we may still accept the representation, if it passes a
</span><span class="comment">*</span><span class="comment">     refined test for RRR. This test supposes that no NaN occurred.
</span><span class="comment">*</span><span class="comment">     Moreover, we use the refined RRR test only for isolated clusters.
</span>      IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND.
     $   (MIN(MAX1,MAX2).LT.FAIL2)
     $  .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
         DORRR1 = .TRUE.
      ELSE
         DORRR1 = .FALSE.
      ENDIF
      TRYRRR1 = .TRUE.
      IF( TRYRRR1 .AND. DORRR1 ) THEN
      IF(INDX.EQ.1) THEN
         TMP = ABS( DPLUS( N ) )
         ZNM2 = ONE
         PROD = ONE
         OLDP = ONE
         DO 15 I = N-1, 1, -1
            IF( PROD .LE. EPS ) THEN
               PROD =
     $         ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
            ELSE
               PROD = PROD*ABS(WORK(N+I))
            END IF
            OLDP = PROD
            ZNM2 = ZNM2 + PROD**2
            TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
 15      CONTINUE
         RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
         IF (RRR1.LE.MAXGROWTH2) THEN
            SIGMA = LSIGMA
            SHIFT = SLEFT
            GOTO 100
         ENDIF
      ELSE IF(INDX.EQ.2) THEN
         TMP = ABS( WORK( N ) )
         ZNM2 = ONE
         PROD = ONE
         OLDP = ONE
         DO 16 I = N-1, 1, -1
            IF( PROD .LE. EPS ) THEN
               PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
            ELSE
               PROD = PROD*ABS(LPLUS(I))
            END IF
            OLDP = PROD
            ZNM2 = ZNM2 + PROD**2
            TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
 16      CONTINUE
         RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
         IF (RRR2.LE.MAXGROWTH2) THEN
            SIGMA = RSIGMA
            SHIFT = SRIGHT
            GOTO 100
         ENDIF
      END IF
      ENDIF

 50   CONTINUE

      IF (KTRY.LT.KTRYMAX) THEN
<span class="comment">*</span><span class="comment">        If we are here, both shifts failed also the RRR test.
</span><span class="comment">*</span><span class="comment">        Back off to the outside
</span>         LSIGMA = MAX( LSIGMA - LDELTA,
     $     LSIGMA - LDMAX)
         RSIGMA = MIN( RSIGMA + RDELTA,
     $     RSIGMA + RDMAX )
         LDELTA = TWO * LDELTA
         RDELTA = TWO * RDELTA
         KTRY = KTRY + 1
         GOTO 5
      ELSE
<span class="comment">*</span><span class="comment">        None of the representations investigated satisfied our
</span><span class="comment">*</span><span class="comment">        criteria. Take the best one we found.
</span>         IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
            LSIGMA = BESTSHIFT
            RSIGMA = BESTSHIFT
            FORCER = .TRUE.
            GOTO 5
         ELSE
            INFO = 1
            RETURN
         ENDIF
      END IF

 100  CONTINUE
      IF (SHIFT.EQ.SLEFT) THEN
      ELSEIF (SHIFT.EQ.SRIGHT) THEN
<span class="comment">*</span><span class="comment">        store new L and D back into DPLUS, LPLUS
</span>         CALL DCOPY( N, WORK, 1, DPLUS, 1 )
         CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
      ENDIF

      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLARRF.371"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>