dstein.f.html

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

</span><span class="comment">*</span><span class="comment">     Get machine constants.
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.177"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize seed for random number generator <a name="DLARNV.179"></a><a href="dlarnv.f.html#DLARNV.1">DLARNV</a>.
</span><span class="comment">*</span><span class="comment">
</span>      DO 40 I = 1, 4
         ISEED( I ) = 1
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize pointers.
</span><span class="comment">*</span><span class="comment">
</span>      INDRV1 = 0
      INDRV2 = INDRV1 + N
      INDRV3 = INDRV2 + N
      INDRV4 = INDRV3 + N
      INDRV5 = INDRV4 + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute eigenvectors of matrix blocks.
</span><span class="comment">*</span><span class="comment">
</span>      J1 = 1
      DO 160 NBLK = 1, IBLOCK( M )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Find starting and ending indices of block nblk.
</span><span class="comment">*</span><span class="comment">
</span>         IF( NBLK.EQ.1 ) THEN
            B1 = 1
         ELSE
            B1 = ISPLIT( NBLK-1 ) + 1
         END IF
         BN = ISPLIT( NBLK )
         BLKSIZ = BN - B1 + 1
         IF( BLKSIZ.EQ.1 )
     $      GO TO 60
         GPIND = B1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute reorthogonalization criterion and stopping criterion.
</span><span class="comment">*</span><span class="comment">
</span>         ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
         ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
         DO 50 I = B1 + 1, BN - 1
            ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+
     $               ABS( E( I ) ) )
   50    CONTINUE
         ORTOL = ODM3*ONENRM
<span class="comment">*</span><span class="comment">
</span>         DTPCRT = SQRT( ODM1 / BLKSIZ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Loop through eigenvalues of block nblk.
</span><span class="comment">*</span><span class="comment">
</span>   60    CONTINUE
         JBLK = 0
         DO 150 J = J1, M
            IF( IBLOCK( J ).NE.NBLK ) THEN
               J1 = J
               GO TO 160
            END IF
            JBLK = JBLK + 1
            XJ = W( J )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Skip all the work if the block size is one.
</span><span class="comment">*</span><span class="comment">
</span>            IF( BLKSIZ.EQ.1 ) THEN
               WORK( INDRV1+1 ) = ONE
               GO TO 120
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If eigenvalues j and j-1 are too close, add a relatively
</span><span class="comment">*</span><span class="comment">           small perturbation.
</span><span class="comment">*</span><span class="comment">
</span>            IF( JBLK.GT.1 ) THEN
               EPS1 = ABS( EPS*XJ )
               PERTOL = TEN*EPS1
               SEP = XJ - XJM
               IF( SEP.LT.PERTOL )
     $            XJ = XJM + PERTOL
            END IF
<span class="comment">*</span><span class="comment">
</span>            ITS = 0
            NRMCHK = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Get random starting vector.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLARNV.258"></a><a href="dlarnv.f.html#DLARNV.1">DLARNV</a>( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy the matrix T so it won't be destroyed in factorization.
</span><span class="comment">*</span><span class="comment">
</span>            CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
            CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
            CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute LU factors with partial pivoting  ( PT = LU )
</span><span class="comment">*</span><span class="comment">
</span>            TOL = ZERO
            CALL <a name="DLAGTF.269"></a><a href="dlagtf.f.html#DLAGTF.1">DLAGTF</a>( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ),
     $                   WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK,
     $                   IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update iteration count.
</span><span class="comment">*</span><span class="comment">
</span>   70       CONTINUE
            ITS = ITS + 1
            IF( ITS.GT.MAXITS )
     $         GO TO 100
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Normalize and scale the righthand side vector Pb.
</span><span class="comment">*</span><span class="comment">
</span>            SCL = BLKSIZ*ONENRM*MAX( EPS,
     $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
     $            DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
            CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve the system LU = Pb.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLAGTS.289"></a><a href="dlagts.f.html#DLAGTS.1">DLAGTS</a>( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ),
     $                   WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK,
     $                   WORK( INDRV1+1 ), TOL, IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reorthogonalize by modified Gram-Schmidt if eigenvalues are
</span><span class="comment">*</span><span class="comment">           close enough.
</span><span class="comment">*</span><span class="comment">
</span>            IF( JBLK.EQ.1 )
     $         GO TO 90
            IF( ABS( XJ-XJM ).GT.ORTOL )
     $         GPIND = J
            IF( GPIND.NE.J ) THEN
               DO 80 I = GPIND, J - 1
                  ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ),
     $                  1 )
                  CALL DAXPY( BLKSIZ, ZTR, Z( B1, I ), 1,
     $                        WORK( INDRV1+1 ), 1 )
   80          CONTINUE
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Check the infinity norm of the iterate.
</span><span class="comment">*</span><span class="comment">
</span>   90       CONTINUE
            JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
            NRM = ABS( WORK( INDRV1+JMAX ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Continue for additional iterations after norm reaches
</span><span class="comment">*</span><span class="comment">           stopping criterion.
</span><span class="comment">*</span><span class="comment">
</span>            IF( NRM.LT.DTPCRT )
     $         GO TO 70
            NRMCHK = NRMCHK + 1
            IF( NRMCHK.LT.EXTRA+1 )
     $         GO TO 70
<span class="comment">*</span><span class="comment">
</span>            GO TO 110
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If stopping criterion was not satisfied, update info and
</span><span class="comment">*</span><span class="comment">           store eigenvector number in array ifail.
</span><span class="comment">*</span><span class="comment">
</span>  100       CONTINUE
            INFO = INFO + 1
            IFAIL( INFO ) = J
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Accept iterate as jth eigenvector.
</span><span class="comment">*</span><span class="comment">
</span>  110       CONTINUE
            SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
            JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
            IF( WORK( INDRV1+JMAX ).LT.ZERO )
     $         SCL = -SCL
            CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
  120       CONTINUE
            DO 130 I = 1, N
               Z( I, J ) = ZERO
  130       CONTINUE
            DO 140 I = 1, BLKSIZ
               Z( B1+I-1, J ) = WORK( INDRV1+I )
  140       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Save the shift to check eigenvalue spacing at next
</span><span class="comment">*</span><span class="comment">           iteration.
</span><span class="comment">*</span><span class="comment">
</span>            XJM = XJ
<span class="comment">*</span><span class="comment">
</span>  150    CONTINUE
  160 CONTINUE
<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="DSTEIN.359"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>