dlasd4.f.html

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

HTML
915 行
字号
1234
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           The following TAU is to approximate
</span><span class="comment">*</span><span class="comment">           SIGMA_n^2 - D( N )*D( N )
</span><span class="comment">*</span><span class="comment">
</span>            IF( A.LT.ZERO ) THEN
               TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
            ELSE
               TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           It can be proved that
</span><span class="comment">*</span><span class="comment">           D(N)^2 &lt; D(N)^2+TAU &lt; SIGMA(N)^2 &lt; D(N)^2+RHO/2
</span><span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        The following ETA is to approximate SIGMA_n - D( N )
</span><span class="comment">*</span><span class="comment">
</span>         ETA = TAU / ( D( N )+SQRT( D( N )*D( N )+TAU ) )
<span class="comment">*</span><span class="comment">
</span>         SIGMA = D( N ) + ETA
         DO 30 J = 1, N
            DELTA( J ) = ( D( J )-D( I ) ) - ETA
            WORK( J ) = D( J ) + D( I ) + ETA
   30    CONTINUE
<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 40 J = 1, II
            TEMP = Z( J ) / ( DELTA( J )*WORK( J ) )
            PSI = PSI + Z( J )*TEMP
            DPSI = DPSI + TEMP*TEMP
            ERRETM = ERRETM + PSI
   40    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 )*WORK( 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">        Test for convergence
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( W ).LE.EPS*ERRETM ) THEN
            GO TO 240
         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
         DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
         DTNSQ = WORK( N )*DELTA( N )
         C = W - DTNSQ1*DPSI - DTNSQ*DPHI
         A = ( DTNSQ+DTNSQ1 )*W - DTNSQ*DTNSQ1*( DPSI+DPHI )
         B = DTNSQ*DTNSQ1*W
         IF( C.LT.ZERO )
     $      C = ABS( C )
         IF( C.EQ.ZERO ) THEN
            ETA = RHO - SIGMA*SIGMA
         ELSE 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 = ETA - DTNSQ
         IF( TEMP.GT.RHO )
     $      ETA = RHO + DTNSQ
<span class="comment">*</span><span class="comment">
</span>         TAU = TAU + ETA
         ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
         DO 50 J = 1, N
            DELTA( J ) = DELTA( J ) - ETA
            WORK( J ) = WORK( J ) + ETA
   50    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         SIGMA = SIGMA + 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 ) / ( WORK( 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 ) / ( WORK( 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
               GO TO 240
            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>            DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
            DTNSQ = WORK( N )*DELTA( N )
            C = W - DTNSQ1*DPSI - DTNSQ*DPHI
            A = ( DTNSQ+DTNSQ1 )*W - DTNSQ1*DTNSQ*( DPSI+DPHI )
            B = DTNSQ1*DTNSQ*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 = ETA - DTNSQ
            IF( TEMP.LE.ZERO )
     $         ETA = ETA / TWO
<span class="comment">*</span><span class="comment">
</span>            TAU = TAU + ETA
            ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
            DO 70 J = 1, N
               DELTA( J ) = DELTA( J ) - ETA
               WORK( J ) = WORK( J ) + ETA
   70       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            SIGMA = SIGMA + 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 ) / ( WORK( 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 ) / ( WORK( 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
         GO TO 240
<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>         DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) )
         DELSQ2 = DELSQ / TWO
         TEMP = DELSQ2 / ( D( I )+SQRT( D( I )*D( I )+DELSQ2 ) )
         DO 100 J = 1, N
            WORK( J ) = D( J ) + D( I ) + TEMP
            DELTA( J ) = ( D( J )-D( I ) ) - TEMP
  100    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         PSI = ZERO
         DO 110 J = 1, I - 1
            PSI = PSI + Z( J )*Z( J ) / ( WORK( 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 ) / ( WORK( J )*DELTA( J ) )
  120    CONTINUE
         C = RHOINV + PSI + PHI
         W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) +
     $       Z( IP1 )*Z( IP1 ) / ( WORK( 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)^2 &lt; the ith sigma^2 &lt; (d(i)^2+d(i+1)^2)/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">
1234

dlasd4.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 dlasd4.f.html 源码文件，采用 HTML 编程语言编写，共 915 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

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