dlaed6.f.html

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

            ELSE
               UBD = TAU
            END IF
            IF( ABS( FINIT ).LE.ABS( TEMP ) )
     $         TAU = ZERO
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     get machine parameters for possible scaling to avoid overflow
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     modified by Sven: parameters SMALL1, SMINV1, SMALL2,
</span><span class="comment">*</span><span class="comment">     SMINV2, EPS are not SAVEd anymore between one call to the
</span><span class="comment">*</span><span class="comment">     others but recomputed at each call
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.172"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )
      BASE = <a name="DLAMCH.173"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Base'</span> )
      SMALL1 = BASE**( INT( LOG( <a name="DLAMCH.174"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'SafMin'</span> ) ) / LOG( BASE ) /
     $         THREE ) )
      SMINV1 = ONE / SMALL1
      SMALL2 = SMALL1*SMALL1
      SMINV2 = SMINV1*SMINV1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine if scaling of inputs necessary to avoid overflow
</span><span class="comment">*</span><span class="comment">     when computing 1/TEMP**3
</span><span class="comment">*</span><span class="comment">
</span>      IF( ORGATI ) THEN
         TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
      ELSE
         TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
      END IF
      SCALE = .FALSE.
      IF( TEMP.LE.SMALL1 ) THEN
         SCALE = .TRUE.
         IF( TEMP.LE.SMALL2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale up by power of radix nearest 1/SAFMIN**(2/3)
</span><span class="comment">*</span><span class="comment">
</span>            SCLFAC = SMINV2
            SCLINV = SMALL2
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale up by power of radix nearest 1/SAFMIN**(1/3)
</span><span class="comment">*</span><span class="comment">
</span>            SCLFAC = SMINV1
            SCLINV = SMALL1
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
</span><span class="comment">*</span><span class="comment">
</span>         DO 10 I = 1, 3
            DSCALE( I ) = D( I )*SCLFAC
            ZSCALE( I ) = Z( I )*SCLFAC
   10    CONTINUE
         TAU = TAU*SCLFAC
         LBD = LBD*SCLFAC
         UBD = UBD*SCLFAC
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy D and Z to DSCALE and ZSCALE
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 I = 1, 3
            DSCALE( I ) = D( I )
            ZSCALE( I ) = Z( I )
   20    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      FC = ZERO
      DF = ZERO
      DDF = ZERO
      DO 30 I = 1, 3
         TEMP = ONE / ( DSCALE( I )-TAU )
         TEMP1 = ZSCALE( I )*TEMP
         TEMP2 = TEMP1*TEMP
         TEMP3 = TEMP2*TEMP
         FC = FC + TEMP1 / DSCALE( I )
         DF = DF + TEMP2
         DDF = DDF + TEMP3
   30 CONTINUE
      F = FINIT + TAU*FC
<span class="comment">*</span><span class="comment">
</span>      IF( ABS( F ).LE.ZERO )
     $   GO TO 60
      IF( F .LE. ZERO )THEN
         LBD = TAU
      ELSE
         UBD = TAU
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
</span><span class="comment">*</span><span class="comment">                            scheme
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     It is not hard to see that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           1) Iterations will go up monotonically
</span><span class="comment">*</span><span class="comment">              if FINIT &lt; 0;
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           2) Iterations will go down monotonically
</span><span class="comment">*</span><span class="comment">              if FINIT &gt; 0.
</span><span class="comment">*</span><span class="comment">
</span>      ITER = NITER + 1
<span class="comment">*</span><span class="comment">
</span>      DO 50 NITER = ITER, MAXIT
<span class="comment">*</span><span class="comment">
</span>         IF( ORGATI ) THEN
            TEMP1 = DSCALE( 2 ) - TAU
            TEMP2 = DSCALE( 3 ) - TAU
         ELSE
            TEMP1 = DSCALE( 1 ) - TAU
            TEMP2 = DSCALE( 2 ) - TAU
         END IF
         A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
         B = TEMP1*TEMP2*F
         C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
         TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
         A = A / TEMP
         B = B / TEMP
         C = C / TEMP
         IF( C.EQ.ZERO ) THEN
            ETA = B / A
         ELSE IF( A.LE.ZERO ) THEN
            ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
         ELSE
            ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
         END IF
         IF( F*ETA.GE.ZERO ) THEN
            ETA = -F / DF
         END IF
<span class="comment">*</span><span class="comment">
</span>         TAU = TAU + ETA
         IF( TAU .LT. LBD .OR. TAU .GT. UBD )
     $      TAU = ( LBD + UBD )/TWO 
<span class="comment">*</span><span class="comment">
</span>         FC = ZERO
         ERRETM = ZERO
         DF = ZERO
         DDF = ZERO
         DO 40 I = 1, 3
            TEMP = ONE / ( DSCALE( I )-TAU )
            TEMP1 = ZSCALE( I )*TEMP
            TEMP2 = TEMP1*TEMP
            TEMP3 = TEMP2*TEMP
            TEMP4 = TEMP1 / DSCALE( I )
            FC = FC + TEMP4
            ERRETM = ERRETM + ABS( TEMP4 )
            DF = DF + TEMP2
            DDF = DDF + TEMP3
   40    CONTINUE
         F = FINIT + TAU*FC
         ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
     $            ABS( TAU )*DF
         IF( ABS( F ).LE.EPS*ERRETM )
     $      GO TO 60
         IF( F .LE. ZERO )THEN
            LBD = TAU
         ELSE
            UBD = TAU
         END IF
   50 CONTINUE
      INFO = 1
   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling
</span><span class="comment">*</span><span class="comment">
</span>      IF( SCALE )
     $   TAU = TAU*SCLINV
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLAED6.325"></a><a href="dlaed6.f.html#DLAED6.1">DLAED6</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>