dtrevc.f.html

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

                           SCALE = SCALE / XNORM
                        END IF
                     END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Scale if necessary
</span><span class="comment">*</span><span class="comment">
</span>                     IF( SCALE.NE.ONE )
     $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
                     WORK( J-1+N ) = X( 1, 1 )
                     WORK( J+N ) = X( 2, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Update right-hand side
</span><span class="comment">*</span><span class="comment">
</span>                     CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
     $                           WORK( 1+N ), 1 )
                     CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
     $                           WORK( 1+N ), 1 )
                  END IF
   60          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy the vector x or Q*x to VR and normalize.
</span><span class="comment">*</span><span class="comment">
</span>               IF( .NOT.OVER ) THEN
                  CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  II = IDAMAX( KI, VR( 1, IS ), 1 )
                  REMAX = ONE / ABS( VR( II, IS ) )
                  CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  DO 70 K = KI + 1, N
                     VR( K, IS ) = ZERO
   70             CONTINUE
               ELSE
                  IF( KI.GT.1 )
     $               CALL DGEMV( <span class="string">'N'</span>, N, KI-1, ONE, VR, LDVR,
     $                           WORK( 1+N ), 1, WORK( KI+N ),
     $                           VR( 1, KI ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  II = IDAMAX( N, VR( 1, KI ), 1 )
                  REMAX = ONE / ABS( VR( II, KI ) )
                  CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
               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">              Complex right eigenvector.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Initial solve
</span><span class="comment">*</span><span class="comment">                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
</span><span class="comment">*</span><span class="comment">                [ (T(KI,KI-1)   T(KI,KI)   )               ]
</span><span class="comment">*</span><span class="comment">
</span>               IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
                  WORK( KI-1+N ) = ONE
                  WORK( KI+N2 ) = WI / T( KI-1, KI )
               ELSE
                  WORK( KI-1+N ) = -WI / T( KI, KI-1 )
                  WORK( KI+N2 ) = ONE
               END IF
               WORK( KI+N ) = ZERO
               WORK( KI-1+N2 ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Form right-hand side
</span><span class="comment">*</span><span class="comment">
</span>               DO 80 K = 1, KI - 2
                  WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
                  WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
   80          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Solve upper quasi-triangular system:
</span><span class="comment">*</span><span class="comment">              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
</span><span class="comment">*</span><span class="comment">
</span>               JNXT = KI - 2
               DO 90 J = KI - 2, 1, -1
                  IF( J.GT.JNXT )
     $               GO TO 90
                  J1 = J
                  J2 = J
                  JNXT = J - 1
                  IF( J.GT.1 ) THEN
                     IF( T( J, J-1 ).NE.ZERO ) THEN
                        J1 = J - 1
                        JNXT = J - 2
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( J1.EQ.J2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    1-by-1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="DLALN2.473"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
     $                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
     $                            X, 2, SCALE, XNORM, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Scale X(1,1) and X(1,2) to avoid overflow when
</span><span class="comment">*</span><span class="comment">                    updating the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                     IF( XNORM.GT.ONE ) THEN
                        IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
                           X( 1, 1 ) = X( 1, 1 ) / XNORM
                           X( 1, 2 ) = X( 1, 2 ) / XNORM
                           SCALE = SCALE / XNORM
                        END IF
                     END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Scale if necessary
</span><span class="comment">*</span><span class="comment">
</span>                     IF( SCALE.NE.ONE ) THEN
                        CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
                        CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
                     END IF
                     WORK( J+N ) = X( 1, 1 )
                     WORK( J+N2 ) = X( 1, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Update the right-hand side
</span><span class="comment">*</span><span class="comment">
</span>                     CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
     $                           WORK( 1+N ), 1 )
                     CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
     $                           WORK( 1+N2 ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    2-by-2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="DLALN2.508"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>( .FALSE., 2, 2, SMIN, ONE,
     $                            T( J-1, J-1 ), LDT, ONE, ONE,
     $                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
     $                            XNORM, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Scale X to avoid overflow when updating
</span><span class="comment">*</span><span class="comment">                    the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                     IF( XNORM.GT.ONE ) THEN
                        BETA = MAX( WORK( J-1 ), WORK( J ) )
                        IF( BETA.GT.BIGNUM / XNORM ) THEN
                           REC = ONE / XNORM
                           X( 1, 1 ) = X( 1, 1 )*REC
                           X( 1, 2 ) = X( 1, 2 )*REC
                           X( 2, 1 ) = X( 2, 1 )*REC
                           X( 2, 2 ) = X( 2, 2 )*REC
                           SCALE = SCALE*REC
                        END IF
                     END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Scale if necessary
</span><span class="comment">*</span><span class="comment">
</span>                     IF( SCALE.NE.ONE ) THEN
                        CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
                        CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
                     END IF
                     WORK( J-1+N ) = X( 1, 1 )
                     WORK( J+N ) = X( 2, 1 )
                     WORK( J-1+N2 ) = X( 1, 2 )
                     WORK( J+N2 ) = X( 2, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Update the right-hand side
</span><span class="comment">*</span><span class="comment">
</span>                     CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
     $                           WORK( 1+N ), 1 )
                     CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
     $                           WORK( 1+N ), 1 )
                     CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
     $                           WORK( 1+N2 ), 1 )
                     CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
     $                           WORK( 1+N2 ), 1 )
                  END IF
   90          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy the vector x or Q*x to VR and normalize.
</span><span class="comment">*</span><span class="comment">
</span>               IF( .NOT.OVER ) THEN
                  CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
                  CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  EMAX = ZERO
                  DO 100 K = 1, KI
                     EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
     $                      ABS( VR( K, IS ) ) )
  100             CONTINUE
<span class="comment">*</span><span class="comment">
</span>                  REMAX = ONE / EMAX
                  CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
                  CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  DO 110 K = KI + 1, N
                     VR( K, IS-1 ) = ZERO
                     VR( K, IS ) = ZERO
  110             CONTINUE
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span>                  IF( KI.GT.2 ) THEN
                     CALL DGEMV( <span class="string">'N'</span>, N, KI-2, ONE, VR, LDVR,
     $                           WORK( 1+N ), 1, WORK( KI-1+N ),
     $                           VR( 1, KI-1 ), 1 )
                     CALL DGEMV( <span class="string">'N'</span>, N, KI-2, ONE, VR, LDVR,
     $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
     $                           VR( 1, KI ), 1 )
                  ELSE
                     CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
                     CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )