ztgevc.f.html

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

                  END IF
                  WORK( J ) = <a name="ZLADIV.423"></a><a href="zladiv.f.html#ZLADIV.1">ZLADIV</a>( -SUM, D )
                  XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
  100          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Back transform eigenvector if HOWMNY='B'.
</span><span class="comment">*</span><span class="comment">
</span>               IF( ILBACK ) THEN
                  CALL ZGEMV( <span class="string">'N'</span>, N, N+1-JE, CONE, VL( 1, JE ), LDVL,
     $                        WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
                  ISRC = 2
                  IBEG = 1
               ELSE
                  ISRC = 1
                  IBEG = JE
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy and scale eigenvector into column of VL
</span><span class="comment">*</span><span class="comment">
</span>               XMAX = ZERO
               DO 110 JR = IBEG, N
                  XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
  110          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( XMAX.GT.SAFMIN ) THEN
                  TEMP = ONE / XMAX
                  DO 120 JR = IBEG, N
                     VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
  120             CONTINUE
               ELSE
                  IBEG = N + 1
               END IF
<span class="comment">*</span><span class="comment">
</span>               DO 130 JR = 1, IBEG - 1
                  VL( JR, IEIG ) = CZERO
  130          CONTINUE
<span class="comment">*</span><span class="comment">
</span>            END IF
  140    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Right eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>      IF( COMPR ) THEN
         IEIG = IM + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Main loop over eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>         DO 250 JE = N, 1, -1
            IF( ILALL ) THEN
               ILCOMP = .TRUE.
            ELSE
               ILCOMP = SELECT( JE )
            END IF
            IF( ILCOMP ) THEN
               IEIG = IEIG - 1
<span class="comment">*</span><span class="comment">
</span>               IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
     $             ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Singular matrix pencil -- return unit eigenvector
</span><span class="comment">*</span><span class="comment">
</span>                  DO 150 JR = 1, N
                     VR( JR, IEIG ) = CZERO
  150             CONTINUE
                  VR( IEIG, IEIG ) = CONE
                  GO TO 250
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Non-singular eigenvalue:
</span><span class="comment">*</span><span class="comment">              Compute coefficients  a  and  b  in
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ( a A - b B ) x  = 0
</span><span class="comment">*</span><span class="comment">
</span>               TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
     $                ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
               SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
               SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
               ACOEFF = SBETA*ASCALE
               BCOEFF = SALPHA*BSCALE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Scale to avoid underflow
</span><span class="comment">*</span><span class="comment">
</span>               LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
               LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
     $               SMALL
<span class="comment">*</span><span class="comment">
</span>               SCALE = ONE
               IF( LSA )
     $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
               IF( LSB )
     $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
     $                    MIN( BNORM, BIG ) )
               IF( LSA .OR. LSB ) THEN
                  SCALE = MIN( SCALE, ONE /
     $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
     $                    ABS1( BCOEFF ) ) ) )
                  IF( LSA ) THEN
                     ACOEFF = ASCALE*( SCALE*SBETA )
                  ELSE
                     ACOEFF = SCALE*ACOEFF
                  END IF
                  IF( LSB ) THEN
                     BCOEFF = BSCALE*( SCALE*SALPHA )
                  ELSE
                     BCOEFF = SCALE*BCOEFF
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               ACOEFA = ABS( ACOEFF )
               BCOEFA = ABS1( BCOEFF )
               XMAX = ONE
               DO 160 JR = 1, N
                  WORK( JR ) = CZERO
  160          CONTINUE
               WORK( JE ) = CONE
               DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Triangular solve of  (a A - b B) x = 0  (columnwise)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              WORK(1:j-1) contains sums w,
</span><span class="comment">*</span><span class="comment">              WORK(j+1:JE) contains x
</span><span class="comment">*</span><span class="comment">
</span>               DO 170 JR = 1, JE - 1
                  WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
  170          CONTINUE
               WORK( JE ) = CONE
<span class="comment">*</span><span class="comment">
</span>               DO 210 J = JE - 1, 1, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Form x(j) := - w(j) / d
</span><span class="comment">*</span><span class="comment">                 with scaling and perturbation of the denominator
</span><span class="comment">*</span><span class="comment">
</span>                  D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
                  IF( ABS1( D ).LE.DMIN )
     $               D = DCMPLX( DMIN )
<span class="comment">*</span><span class="comment">
</span>                  IF( ABS1( D ).LT.ONE ) THEN
                     IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
                        TEMP = ONE / ABS1( WORK( J ) )
                        DO 180 JR = 1, JE
                           WORK( JR ) = TEMP*WORK( JR )
  180                   CONTINUE
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  WORK( J ) = <a name="ZLADIV.568"></a><a href="zladiv.f.html#ZLADIV.1">ZLADIV</a>( -WORK( J ), D )
<span class="comment">*</span><span class="comment">
</span>                  IF( J.GT.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
</span><span class="comment">*</span><span class="comment">
</span>                     IF( ABS1( WORK( J ) ).GT.ONE ) THEN
                        TEMP = ONE / ABS1( WORK( J ) )
                        IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
     $                      BIGNUM*TEMP ) THEN
                           DO 190 JR = 1, JE
                              WORK( JR ) = TEMP*WORK( JR )
  190                      CONTINUE
                        END IF
                     END IF
<span class="comment">*</span><span class="comment">
</span>                     CA = ACOEFF*WORK( J )
                     CB = BCOEFF*WORK( J )
                     DO 200 JR = 1, J - 1
                        WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
     $                               CB*P( JR, J )
  200                CONTINUE
                  END IF
  210          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Back transform eigenvector if HOWMNY='B'.
</span><span class="comment">*</span><span class="comment">
</span>               IF( ILBACK ) THEN
                  CALL ZGEMV( <span class="string">'N'</span>, N, JE, CONE, VR, LDVR, WORK, 1,
     $                        CZERO, WORK( N+1 ), 1 )
                  ISRC = 2
                  IEND = N
               ELSE
                  ISRC = 1
                  IEND = JE
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy and scale eigenvector into column of VR
</span><span class="comment">*</span><span class="comment">
</span>               XMAX = ZERO
               DO 220 JR = 1, IEND
                  XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
  220          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( XMAX.GT.SAFMIN ) THEN
                  TEMP = ONE / XMAX
                  DO 230 JR = 1, IEND
                     VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
  230             CONTINUE
               ELSE
                  IEND = 0
               END IF
<span class="comment">*</span><span class="comment">
</span>               DO 240 JR = IEND + 1, N
                  VR( JR, IEIG ) = CZERO
  240          CONTINUE
<span class="comment">*</span><span class="comment">
</span>            END IF
  250    CONTINUE
      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="ZTGEVC.631"></a><a href="ztgevc.f.html#ZTGEVC.1">ZTGEVC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>