slaebz.f.html

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

                  END IF
   70          CONTINUE
               IF( INFO.NE.0 )
     $            RETURN
               KL = KLNEW
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              IJOB=3: Binary search.  Keep only the interval containing
</span><span class="comment">*</span><span class="comment">                      w   s.t. N(w) = NVAL
</span><span class="comment">*</span><span class="comment">
</span>               DO 80 JI = KF, KL
                  IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
                     AB( JI, 1 ) = C( JI )
                     NAB( JI, 1 ) = IWORK( JI )
                  END IF
                  IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
                     AB( JI, 2 ) = C( JI )
                     NAB( JI, 2 ) = IWORK( JI )
                  END IF
   80          CONTINUE
            END IF
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           End of Parallel Version of the loop
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Begin of Serial Version of the loop
</span><span class="comment">*</span><span class="comment">
</span>            KLNEW = KL
            DO 100 JI = KF, KL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute N(w), the number of eigenvalues less than w
</span><span class="comment">*</span><span class="comment">
</span>               TMP1 = C( JI )
               TMP2 = D( 1 ) - TMP1
               ITMP1 = 0
               IF( TMP2.LE.PIVMIN ) THEN
                  ITMP1 = 1
                  TMP2 = MIN( TMP2, -PIVMIN )
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A series of compiler directives to defeat vectorization
</span><span class="comment">*</span><span class="comment">              for the next loop
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">$PL$ CMCHAR=' '
</span>CDIR$          NEXTSCALAR
C$DIR          SCALAR
CDIR$          NEXT SCALAR
CVD$L          NOVECTOR
CDEC$          NOVECTOR
CVD$           NOVECTOR
<span class="comment">*</span><span class="comment">VDIR          NOVECTOR
</span><span class="comment">*</span><span class="comment">VOCL          LOOP,SCALAR
</span>CIBM           PREFER SCALAR
<span class="comment">*</span><span class="comment">$PL$ CMCHAR='*'
</span><span class="comment">*</span><span class="comment">
</span>               DO 90 J = 2, N
                  TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
                  IF( TMP2.LE.PIVMIN ) THEN
                     ITMP1 = ITMP1 + 1
                     TMP2 = MIN( TMP2, -PIVMIN )
                  END IF
   90          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( IJOB.LE.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 IJOB=2: Choose all intervals containing eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Insure that N(w) is monotone
</span><span class="comment">*</span><span class="comment">
</span>                  ITMP1 = MIN( NAB( JI, 2 ),
     $                    MAX( NAB( JI, 1 ), ITMP1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Update the Queue -- add intervals if both halves
</span><span class="comment">*</span><span class="comment">                 contain eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    No eigenvalue in the upper interval:
</span><span class="comment">*</span><span class="comment">                    just use the lower interval.
</span><span class="comment">*</span><span class="comment">
</span>                     AB( JI, 2 ) = TMP1
<span class="comment">*</span><span class="comment">
</span>                  ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    No eigenvalue in the lower interval:
</span><span class="comment">*</span><span class="comment">                    just use the upper interval.
</span><span class="comment">*</span><span class="comment">
</span>                     AB( JI, 1 ) = TMP1
                  ELSE IF( KLNEW.LT.MMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Eigenvalue in both intervals -- add upper to queue.
</span><span class="comment">*</span><span class="comment">
</span>                     KLNEW = KLNEW + 1
                     AB( KLNEW, 2 ) = AB( JI, 2 )
                     NAB( KLNEW, 2 ) = NAB( JI, 2 )
                     AB( KLNEW, 1 ) = TMP1
                     NAB( KLNEW, 1 ) = ITMP1
                     AB( JI, 2 ) = TMP1
                     NAB( JI, 2 ) = ITMP1
                  ELSE
                     INFO = MMAX + 1
                     RETURN
                  END IF
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 IJOB=3: Binary search.  Keep only the interval
</span><span class="comment">*</span><span class="comment">                         containing  w  s.t. N(w) = NVAL
</span><span class="comment">*</span><span class="comment">
</span>                  IF( ITMP1.LE.NVAL( JI ) ) THEN
                     AB( JI, 1 ) = TMP1
                     NAB( JI, 1 ) = ITMP1
                  END IF
                  IF( ITMP1.GE.NVAL( JI ) ) THEN
                     AB( JI, 2 ) = TMP1
                     NAB( JI, 2 ) = ITMP1
                  END IF
               END IF
  100       CONTINUE
            KL = KLNEW
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           End of Serial Version of the loop
</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">        Check for convergence
</span><span class="comment">*</span><span class="comment">
</span>         KFNEW = KF
         DO 110 JI = KF, KL
            TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
            TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
            IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
     $          NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Converged -- Swap with position KFNEW,
</span><span class="comment">*</span><span class="comment">                           then increment KFNEW
</span><span class="comment">*</span><span class="comment">
</span>               IF( JI.GT.KFNEW ) THEN
                  TMP1 = AB( JI, 1 )
                  TMP2 = AB( JI, 2 )
                  ITMP1 = NAB( JI, 1 )
                  ITMP2 = NAB( JI, 2 )
                  AB( JI, 1 ) = AB( KFNEW, 1 )
                  AB( JI, 2 ) = AB( KFNEW, 2 )
                  NAB( JI, 1 ) = NAB( KFNEW, 1 )
                  NAB( JI, 2 ) = NAB( KFNEW, 2 )
                  AB( KFNEW, 1 ) = TMP1
                  AB( KFNEW, 2 ) = TMP2
                  NAB( KFNEW, 1 ) = ITMP1
                  NAB( KFNEW, 2 ) = ITMP2
                  IF( IJOB.EQ.3 ) THEN
                     ITMP1 = NVAL( JI )
                     NVAL( JI ) = NVAL( KFNEW )
                     NVAL( KFNEW ) = ITMP1
                  END IF
               END IF
               KFNEW = KFNEW + 1
            END IF
  110    CONTINUE
         KF = KFNEW
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Choose Midpoints
</span><span class="comment">*</span><span class="comment">
</span>         DO 120 JI = KF, KL
            C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If no more intervals to refine, quit.
</span><span class="comment">*</span><span class="comment">
</span>         IF( KF.GT.KL )
     $      GO TO 140
  130 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Converged
</span><span class="comment">*</span><span class="comment">
</span>  140 CONTINUE
      INFO = MAX( KL+1-KF, 0 )
      MOUT = KL
<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="SLAEBZ.549"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>