slaed4.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 869 行 · 第 1/3 页
HTML
869 行
</span><span class="comment">*</span><span class="comment">        if for some reason caused by roundoff, eta*w &gt; 0,
</span><span class="comment">*</span><span class="comment">        we simply use one Newton step instead. This way
</span><span class="comment">*</span><span class="comment">        will guarantee eta*w &lt; 0.
</span><span class="comment">*</span><span class="comment">
</span>         IF( W*ETA.GT.ZERO )
     $      ETA = -W / ( DPSI+DPHI )
         TEMP = TAU + ETA
         IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
            IF( W.LT.ZERO ) THEN
               ETA = ( DLTUB-TAU ) / TWO
            ELSE
               ETA = ( DLTLB-TAU ) / TWO
            END IF
         END IF
         DO 50 J = 1, N
            DELTA( J ) = DELTA( J ) - ETA
   50    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         TAU = TAU + ETA
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Evaluate PSI and the derivative DPSI
</span><span class="comment">*</span><span class="comment">
</span>         DPSI = ZERO
         PSI = ZERO
         ERRETM = ZERO
         DO 60 J = 1, II
            TEMP = Z( J ) / DELTA( J )
            PSI = PSI + Z( J )*TEMP
            DPSI = DPSI + TEMP*TEMP
            ERRETM = ERRETM + PSI
   60    CONTINUE
         ERRETM = ABS( ERRETM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Evaluate PHI and the derivative DPHI
</span><span class="comment">*</span><span class="comment">
</span>         TEMP = Z( N ) / DELTA( N )
         PHI = Z( N )*TEMP
         DPHI = TEMP*TEMP
         ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
     $            ABS( TAU )*( DPSI+DPHI )
<span class="comment">*</span><span class="comment">
</span>         W = RHOINV + PHI + PSI
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Main loop to update the values of the array   DELTA
</span><span class="comment">*</span><span class="comment">
</span>         ITER = NITER + 1
<span class="comment">*</span><span class="comment">
</span>         DO 90 NITER = ITER, MAXIT
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test for convergence
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( W ).LE.EPS*ERRETM ) THEN
               DLAM = D( I ) + TAU
               GO TO 250
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( W.LE.ZERO ) THEN
               DLTLB = MAX( DLTLB, TAU )
            ELSE
               DLTUB = MIN( DLTUB, TAU )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Calculate the new step
</span><span class="comment">*</span><span class="comment">
</span>            C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
            A = ( DELTA( N-1 )+DELTA( N ) )*W -
     $          DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
            B = DELTA( N-1 )*DELTA( N )*W
            IF( A.GE.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
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Note, eta should be positive if w is negative, and
</span><span class="comment">*</span><span class="comment">           eta should be negative otherwise. However,
</span><span class="comment">*</span><span class="comment">           if for some reason caused by roundoff, eta*w &gt; 0,
</span><span class="comment">*</span><span class="comment">           we simply use one Newton step instead. This way
</span><span class="comment">*</span><span class="comment">           will guarantee eta*w &lt; 0.
</span><span class="comment">*</span><span class="comment">
</span>            IF( W*ETA.GT.ZERO )
     $         ETA = -W / ( DPSI+DPHI )
            TEMP = TAU + ETA
            IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
               IF( W.LT.ZERO ) THEN
                  ETA = ( DLTUB-TAU ) / TWO
               ELSE
                  ETA = ( DLTLB-TAU ) / TWO
               END IF
            END IF
            DO 70 J = 1, N
               DELTA( J ) = DELTA( J ) - ETA
   70       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            TAU = TAU + ETA
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Evaluate PSI and the derivative DPSI
</span><span class="comment">*</span><span class="comment">
</span>            DPSI = ZERO
            PSI = ZERO
            ERRETM = ZERO
            DO 80 J = 1, II
               TEMP = Z( J ) / DELTA( J )
               PSI = PSI + Z( J )*TEMP
               DPSI = DPSI + TEMP*TEMP
               ERRETM = ERRETM + PSI
   80       CONTINUE
            ERRETM = ABS( ERRETM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Evaluate PHI and the derivative DPHI
</span><span class="comment">*</span><span class="comment">
</span>            TEMP = Z( N ) / DELTA( N )
            PHI = Z( N )*TEMP
            DPHI = TEMP*TEMP
            ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
     $               ABS( TAU )*( DPSI+DPHI )
<span class="comment">*</span><span class="comment">
</span>            W = RHOINV + PHI + PSI
   90    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Return with INFO = 1, NITER = MAXIT and not converged
</span><span class="comment">*</span><span class="comment">
</span>         INFO = 1
         DLAM = D( I ) + TAU
         GO TO 250
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        End for the case I = N
</span><span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        The case for I &lt; N
</span><span class="comment">*</span><span class="comment">
</span>         NITER = 1
         IP1 = I + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Calculate initial guess
</span><span class="comment">*</span><span class="comment">
</span>         DEL = D( IP1 ) - D( I )
         MIDPT = DEL / TWO
         DO 100 J = 1, N
            DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
  100    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         PSI = ZERO
         DO 110 J = 1, I - 1
            PSI = PSI + Z( J )*Z( J ) / DELTA( J )
  110    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         PHI = ZERO
         DO 120 J = N, I + 2, -1
            PHI = PHI + Z( J )*Z( J ) / DELTA( J )
  120    CONTINUE
         C = RHOINV + PSI + PHI
         W = C + Z( I )*Z( I ) / DELTA( I ) +
     $       Z( IP1 )*Z( IP1 ) / DELTA( IP1 )
<span class="comment">*</span><span class="comment">
</span>         IF( W.GT.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           d(i)&lt; the ith eigenvalue &lt; (d(i)+d(i+1))/2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           We choose d(i) as origin.
</span><span class="comment">*</span><span class="comment">
</span>            ORGATI = .TRUE.
            A = C*DEL + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
            B = Z( I )*Z( I )*DEL
            IF( A.GT.ZERO ) THEN
               TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
            ELSE
               TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
            END IF
            DLTLB = ZERO
            DLTUB = MIDPT
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           (d(i)+d(i+1))/2 &lt;= the ith eigenvalue &lt; d(i+1)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           We choose d(i+1) as origin.
</span><span class="comment">*</span><span class="comment">
</span>            ORGATI = .FALSE.
            A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
            B = Z( IP1 )*Z( IP1 )*DEL
            IF( A.LT.ZERO ) THEN
               TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
            ELSE
               TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
            END IF
            DLTLB = -MIDPT
            DLTUB = ZERO
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( ORGATI ) THEN
            DO 130 J = 1, N
               DELTA( J ) = ( D( J )-D( I ) ) - TAU
  130       CONTINUE
         ELSE
            DO 140 J = 1, N
               DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU
  140       CONTINUE
         END IF
         IF( ORGATI ) THEN
            II = I
         ELSE
            II = I + 1
         END IF
         IIM1 = II - 1
         IIP1 = II + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Evaluate PSI and the derivative DPSI
</span><span class="comment">*</span><span class="comment">
</span>         DPSI = ZERO
         PSI = ZERO
         ERRETM = ZERO
         DO 150 J = 1, IIM1
            TEMP = Z( J ) / DELTA( J )
            PSI = PSI + Z( J )*TEMP
            DPSI = DPSI + TEMP*TEMP
            ERRETM = ERRETM + PSI
  150    CONTINUE
         ERRETM = ABS( ERRETM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Evaluate PHI and the derivative DPHI
</span><span class="comment">*</span><span class="comment">
</span>         DPHI = ZERO
         PHI = ZERO
         DO 160 J = N, IIP1, -1
            TEMP = Z( J ) / DELTA( J )
            PHI = PHI + Z( J )*TEMP
            DPHI = DPHI + TEMP*TEMP
            ERRETM = ERRETM + PHI
  160    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         W = RHOINV + PHI + PSI
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        W is the value of the secular function with
</span><span class="comment">*</span><span class="comment">        its ii-th element removed.
</span><span class="comment">*</span><span class="comment">
</span>         SWTCH3 = .FALSE.
         IF( ORGATI ) THEN
            IF( W.LT.ZERO )
     $         SWTCH3 = .TRUE.
         ELSE
            IF( W.GT.ZERO )
     $         SWTCH3 = .TRUE.
         END IF
         IF( II.EQ.1 .OR. II.EQ.N )
     $      SWTCH3 = .FALSE.
<span class="comment">*</span><span class="comment">
</span>         TEMP = Z( II ) / DELTA( II )
         DW = DPSI + DPHI + TEMP*TEMP
         TEMP = Z( II )*TEMP
         W = W + TEMP
         ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
     $            THREE*ABS( TEMP ) + ABS( TAU )*DW
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Test for convergence
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( W ).LE.EPS*ERRETM ) THEN
            IF( ORGATI ) THEN
               DLAM = D( I ) + TAU
            ELSE
               DLAM = D( IP1 ) + TAU
            END IF
            GO TO 250
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( W.LE.ZERO ) THEN
            DLTLB = MAX( DLTLB, TAU )
         ELSE
            DLTUB = MIN( DLTUB, TAU )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Calculate the new step
</span><span class="comment">*</span><span class="comment">
</span>         NITER = NITER + 1
         IF( .NOT.SWTCH3 ) THEN
            IF( ORGATI ) THEN
               C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )*
     $             ( Z( I ) / DELTA( I ) )**2
            ELSE
               C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )*
     $             ( Z( IP1 ) / DELTA( IP1 ) )**2
            END IF
            A = ( DELTA( I )+DELTA( IP1 ) )*W -
     $          DELTA( I )*DELTA( IP1 )*DW
            B = DELTA( I )*DELTA( IP1 )*W
            IF( C.EQ.ZERO ) THEN
               IF( A.EQ.ZERO ) THEN
                  IF( ORGATI ) THEN
                     A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )*
     $                   ( DPSI+DPHI )
slaed4.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 slaed4.f.html 源码文件，采用 HTML 编程语言编写，共 869 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。
虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。
⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?