slaein.f.html

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

            I3 = -1
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           UL decomposition with partial pivoting of conjg(B),
</span><span class="comment">*</span><span class="comment">           replacing zero 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( N+1, N ) = WI
            DO 180 J = 1, N - 1
               B( N+1, J ) = ZERO
  180       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            DO 210 J = N, 2, -1
               EJ = H( J, J-1 )
               ABSBJJ = <a name="SLAPY2.370"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( B( J, J ), B( J+1, J ) )
               IF( ABSBJJ.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>                  XR = B( J, J ) / EJ
                  XI = B( J+1, J ) / EJ
                  B( J, J ) = EJ
                  B( J+1, J ) = ZERO
                  DO 190 I = 1, J - 1
                     TEMP = B( I, J-1 )
                     B( I, J-1 ) = B( I, J ) - XR*TEMP
                     B( J, I ) = B( J+1, I ) - XI*TEMP
                     B( I, J ) = TEMP
                     B( J+1, I ) = ZERO
  190             CONTINUE
                  B( J+1, J-1 ) = WI
                  B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
                  B( J, J-1 ) = B( J, J-1 ) - XR*WI
               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( ABSBJJ.EQ.ZERO ) THEN
                     B( J, J ) = EPS3
                     B( J+1, J ) = ZERO
                     ABSBJJ = EPS3
                  END IF
                  EJ = ( EJ / ABSBJJ ) / ABSBJJ
                  XR = B( J, J )*EJ
                  XI = -B( J+1, J )*EJ
                  DO 200 I = 1, J - 1
                     B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
     $                             XI*B( J+1, I )
                     B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
  200             CONTINUE
                  B( J, J-1 ) = B( J, J-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 j-th column.
</span><span class="comment">*</span><span class="comment">
</span>               WORK( J ) = SASUM( J-1, B( 1, J ), 1 ) +
     $                     SASUM( J-1, B( J+1, 1 ), LDB )
  210       CONTINUE
            IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
     $         B( 1, 1 ) = EPS3
            WORK( 1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>            I1 = 1
            I2 = N
            I3 = 1
         END IF
<span class="comment">*</span><span class="comment">
</span>         DO 270 ITS = 1, N
            SCALE = ONE
            VMAX = ONE
            VCRIT = BIGNUM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
</span><span class="comment">*</span><span class="comment">             or U'*(xr,xi) = scale*(vr,vi) for a left eigenvector,
</span><span class="comment">*</span><span class="comment">           overwriting (xr,xi) on (vr,vi).
</span><span class="comment">*</span><span class="comment">
</span>            DO 250 I = I1, I2, I3
<span class="comment">*</span><span class="comment">
</span>               IF( WORK( I ).GT.VCRIT ) THEN
                  REC = ONE / VMAX
                  CALL SSCAL( N, REC, VR, 1 )
                  CALL SSCAL( N, REC, VI, 1 )
                  SCALE = SCALE*REC
                  VMAX = ONE
                  VCRIT = BIGNUM
               END IF
<span class="comment">*</span><span class="comment">
</span>               XR = VR( I )
               XI = VI( I )
               IF( RIGHTV ) THEN
                  DO 220 J = I + 1, N
                     XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
                     XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
  220             CONTINUE
               ELSE
                  DO 230 J = 1, I - 1
                     XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
                     XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
  230             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span>               W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
               IF( W.GT.SMLNUM ) THEN
                  IF( W.LT.ONE ) THEN
                     W1 = ABS( XR ) + ABS( XI )
                     IF( W1.GT.W*BIGNUM ) THEN
                        REC = ONE / W1
                        CALL SSCAL( N, REC, VR, 1 )
                        CALL SSCAL( N, REC, VI, 1 )
                        XR = VR( I )
                        XI = VI( I )
                        SCALE = SCALE*REC
                        VMAX = VMAX*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Divide by diagonal element of B.
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SLADIV.474"></a><a href="sladiv.f.html#SLADIV.1">SLADIV</a>( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
     $                         VI( I ) )
                  VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
                  VCRIT = BIGNUM / VMAX
               ELSE
                  DO 240 J = 1, N
                     VR( J ) = ZERO
                     VI( J ) = ZERO
  240             CONTINUE
                  VR( I ) = ONE
                  VI( I ) = ONE
                  SCALE = ZERO
                  VMAX = ONE
                  VCRIT = BIGNUM
               END IF
  250       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test for sufficient growth in the norm of (VR,VI).
</span><span class="comment">*</span><span class="comment">
</span>            VNORM = SASUM( N, VR, 1 ) + SASUM( N, VI, 1 )
            IF( VNORM.GE.GROWTO*SCALE )
     $         GO TO 280
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Choose a new orthogonal starting vector and try again.
</span><span class="comment">*</span><span class="comment">
</span>            Y = EPS3 / ( ROOTN+ONE )
            VR( 1 ) = EPS3
            VI( 1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>            DO 260 I = 2, N
               VR( I ) = Y
               VI( I ) = ZERO
  260       CONTINUE
            VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
  270    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>  280    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>         VNORM = ZERO
         DO 290 I = 1, N
            VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
  290    CONTINUE
         CALL SSCAL( N, ONE / VNORM, VR, 1 )
         CALL SSCAL( N, ONE / VNORM, VI, 1 )
<span class="comment">*</span><span class="comment">
</span>      END IF
<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="SLAEIN.529"></a><a href="slaein.f.html#SLAEIN.1">SLAEIN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>