slaein.f.html

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

                     TEMP = B( I+1, J )
                     B( I+1, J ) = B( I, J ) - X*TEMP
                     B( I, J ) = TEMP
   40             CONTINUE
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Eliminate without interchange.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( B( I, I ).EQ.ZERO )
     $               B( I, I ) = EPS3
                  X = EI / B( I, I )
                  IF( X.NE.ZERO ) THEN
                     DO 50 J = I + 1, N
                        B( I+1, J ) = B( I+1, J ) - X*B( I, J )
   50                CONTINUE
                  END IF
               END IF
   60       CONTINUE
            IF( B( N, N ).EQ.ZERO )
     $         B( N, N ) = EPS3
<span class="comment">*</span><span class="comment">
</span>            TRANS = <span class="string">'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">           UL decomposition with partial pivoting of B, replacing zero
</span><span class="comment">*</span><span class="comment">           pivots by EPS3.
</span><span class="comment">*</span><span class="comment">
</span>            DO 90 J = N, 2, -1
               EJ = H( J, J-1 )
               IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Interchange columns and eliminate.
</span><span class="comment">*</span><span class="comment">
</span>                  X = B( J, J ) / EJ
                  B( J, J ) = EJ
                  DO 70 I = 1, J - 1
                     TEMP = B( I, J-1 )
                     B( I, J-1 ) = B( I, J ) - X*TEMP
                     B( I, J ) = TEMP
   70             CONTINUE
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Eliminate without interchange.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( B( J, J ).EQ.ZERO )
     $               B( J, J ) = EPS3
                  X = EJ / B( J, J )
                  IF( X.NE.ZERO ) THEN
                     DO 80 I = 1, J - 1
                        B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
   80                CONTINUE
                  END IF
               END IF
   90       CONTINUE
            IF( B( 1, 1 ).EQ.ZERO )
     $         B( 1, 1 ) = EPS3
<span class="comment">*</span><span class="comment">
</span>            TRANS = <span class="string">'T'</span>
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span>         NORMIN = <span class="string">'N'</span>
         DO 110 ITS = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve U*x = scale*v for a right eigenvector
</span><span class="comment">*</span><span class="comment">             or U'*x = scale*v for a left eigenvector,
</span><span class="comment">*</span><span class="comment">           overwriting x on v.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLATRS.237"></a><a href="slatrs.f.html#SLATRS.1">SLATRS</a>( <span class="string">'Upper'</span>, TRANS, <span class="string">'Nonunit'</span>, NORMIN, N, B, LDB,
     $                   VR, SCALE, WORK, IERR )
            NORMIN = <span class="string">'Y'</span>
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test for sufficient growth in the norm of v.
</span><span class="comment">*</span><span class="comment">
</span>            VNORM = SASUM( N, VR, 1 )
            IF( VNORM.GE.GROWTO*SCALE )
     $         GO TO 120
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Choose new orthogonal starting vector and try again.
</span><span class="comment">*</span><span class="comment">
</span>            TEMP = EPS3 / ( ROOTN+ONE )
            VR( 1 ) = EPS3
            DO 100 I = 2, N
               VR( I ) = TEMP
  100       CONTINUE
            VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
  110    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Failure to find eigenvector in N iterations.
</span><span class="comment">*</span><span class="comment">
</span>         INFO = 1
<span class="comment">*</span><span class="comment">
</span>  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Normalize eigenvector.
</span><span class="comment">*</span><span class="comment">
</span>         I = ISAMAX( N, VR, 1 )
         CALL SSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Complex eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>         IF( NOINIT ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Set initial vector.
</span><span class="comment">*</span><span class="comment">
</span>            DO 130 I = 1, N
               VR( I ) = EPS3
               VI( I ) = ZERO
  130       CONTINUE
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Scale supplied initial vector.
</span><span class="comment">*</span><span class="comment">
</span>            NORM = <a name="SLAPY2.283"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( SNRM2( N, VR, 1 ), SNRM2( N, VI, 1 ) )
            REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
            CALL SSCAL( N, REC, VR, 1 )
            CALL SSCAL( N, REC, VI, 1 )
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( RIGHTV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           LU decomposition with partial pivoting of B, replacing zero
</span><span class="comment">*</span><span class="comment">           pivots by EPS3.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           The imaginary part of the (i,j)-th element of U is stored in
</span><span class="comment">*</span><span class="comment">           B(j+1,i).
</span><span class="comment">*</span><span class="comment">
</span>            B( 2, 1 ) = -WI
            DO 140 I = 2, N
               B( I+1, 1 ) = ZERO
  140       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            DO 170 I = 1, N - 1
               ABSBII = <a name="SLAPY2.303"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( B( I, I ), B( I+1, I ) )
               EI = H( I+1, I )
               IF( ABSBII.LT.ABS( EI ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Interchange rows and eliminate.
</span><span class="comment">*</span><span class="comment">
</span>                  XR = B( I, I ) / EI
                  XI = B( I+1, I ) / EI
                  B( I, I ) = EI
                  B( I+1, I ) = ZERO
                  DO 150 J = I + 1, N
                     TEMP = B( I+1, J )
                     B( I+1, J ) = B( I, J ) - XR*TEMP
                     B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
                     B( I, J ) = TEMP
                     B( J+1, I ) = ZERO
  150             CONTINUE
                  B( I+2, I ) = -WI
                  B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
                  B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Eliminate without interchanging rows.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( ABSBII.EQ.ZERO ) THEN
                     B( I, I ) = EPS3
                     B( I+1, I ) = ZERO
                     ABSBII = EPS3
                  END IF
                  EI = ( EI / ABSBII ) / ABSBII
                  XR = B( I, I )*EI
                  XI = -B( I+1, I )*EI
                  DO 160 J = I + 1, N
                     B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
     $                             XI*B( J+1, I )
                     B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
  160             CONTINUE
                  B( I+2, I+1 ) = B( I+2, I+1 ) - WI
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute 1-norm of offdiagonal elements of i-th row.
</span><span class="comment">*</span><span class="comment">
</span>               WORK( I ) = SASUM( N-I, B( I, I+1 ), LDB ) +
     $                     SASUM( N-I, B( I+2, I ), 1 )
  170       CONTINUE
            IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
     $         B( N, N ) = EPS3
            WORK( N ) = ZERO
<span class="comment">*</span><span class="comment">
</span>            I1 = N
            I2 = 1