dtrevc.f.html

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

               ELSE
                  WORK( KI+N ) = ONE
                  WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
               END IF
               WORK( KI+1+N ) = ZERO
               WORK( KI+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 190 K = KI + 2, N
                  WORK( K+N ) = -WORK( KI+N )*T( KI, K )
                  WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
  190          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Solve complex quasi-triangular system:
</span><span class="comment">*</span><span class="comment">              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
</span><span class="comment">*</span><span class="comment">
</span>               VMAX = ONE
               VCRIT = BIGNUM
<span class="comment">*</span><span class="comment">
</span>               JNXT = KI + 2
               DO 200 J = KI + 2, N
                  IF( J.LT.JNXT )
     $               GO TO 200
                  J1 = J
                  J2 = J
                  JNXT = J + 1
                  IF( J.LT.N ) THEN
                     IF( T( J+1, J ).NE.ZERO ) THEN
                        J2 = 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><span class="comment">*</span><span class="comment">                    Scale if necessary to avoid overflow when
</span><span class="comment">*</span><span class="comment">                    forming the right-hand side elements.
</span><span class="comment">*</span><span class="comment">
</span>                     IF( WORK( J ).GT.VCRIT ) THEN
                        REC = ONE / VMAX
                        CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
                        CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
                        VMAX = ONE
                        VCRIT = BIGNUM
                     END IF
<span class="comment">*</span><span class="comment">
</span>                     WORK( J+N ) = WORK( J+N ) -
     $                             DDOT( J-KI-2, T( KI+2, J ), 1,
     $                             WORK( KI+2+N ), 1 )
                     WORK( J+N2 ) = WORK( J+N2 ) -
     $                              DDOT( J-KI-2, T( KI+2, J ), 1,
     $                              WORK( KI+2+N2 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="DLALN2.844"></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 if necessary
</span><span class="comment">*</span><span class="comment">
</span>                     IF( SCALE.NE.ONE ) THEN
                        CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
                        CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
                     END IF
                     WORK( J+N ) = X( 1, 1 )
                     WORK( J+N2 ) = X( 1, 2 )
                     VMAX = MAX( ABS( WORK( J+N ) ),
     $                      ABS( WORK( J+N2 ) ), VMAX )
                     VCRIT = BIGNUM / VMAX
<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><span class="comment">*</span><span class="comment">                    Scale if necessary to avoid overflow when forming
</span><span class="comment">*</span><span class="comment">                    the right-hand side elements.
</span><span class="comment">*</span><span class="comment">
</span>                     BETA = MAX( WORK( J ), WORK( J+1 ) )
                     IF( BETA.GT.VCRIT ) THEN
                        REC = ONE / VMAX
                        CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
                        CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
                        VMAX = ONE
                        VCRIT = BIGNUM
                     END IF
<span class="comment">*</span><span class="comment">
</span>                     WORK( J+N ) = WORK( J+N ) -
     $                             DDOT( J-KI-2, T( KI+2, J ), 1,
     $                             WORK( KI+2+N ), 1 )
<span class="comment">*</span><span class="comment">
</span>                     WORK( J+N2 ) = WORK( J+N2 ) -
     $                              DDOT( J-KI-2, T( KI+2, J ), 1,
     $                              WORK( KI+2+N2 ), 1 )
<span class="comment">*</span><span class="comment">
</span>                     WORK( J+1+N ) = WORK( J+1+N ) -
     $                               DDOT( J-KI-2, T( KI+2, J+1 ), 1,
     $                               WORK( KI+2+N ), 1 )
<span class="comment">*</span><span class="comment">
</span>                     WORK( J+1+N2 ) = WORK( J+1+N2 ) -
     $                                DDOT( J-KI-2, T( KI+2, J+1 ), 1,
     $                                WORK( KI+2+N2 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Solve 2-by-2 complex linear equation
</span><span class="comment">*</span><span class="comment">                      ([T(j,j)   T(j,j+1)  ]'-(wr-i*wi)*I)*X = SCALE*B
</span><span class="comment">*</span><span class="comment">                      ([T(j+1,j) T(j+1,j+1)]             )
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="DLALN2.896"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>( .TRUE., 2, 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 if necessary
</span><span class="comment">*</span><span class="comment">
</span>                     IF( SCALE.NE.ONE ) THEN
                        CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
                        CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
                     END IF
                     WORK( J+N ) = X( 1, 1 )
                     WORK( J+N2 ) = X( 1, 2 )
                     WORK( J+1+N ) = X( 2, 1 )
                     WORK( J+1+N2 ) = X( 2, 2 )
                     VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
     $                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
                     VCRIT = BIGNUM / VMAX
<span class="comment">*</span><span class="comment">
</span>                  END IF
  200          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy the vector x or Q*x to VL and normalize.
</span><span class="comment">*</span><span class="comment">
</span>               IF( .NOT.OVER ) THEN
                  CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
                  CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
     $                        1 )
<span class="comment">*</span><span class="comment">
</span>                  EMAX = ZERO
                  DO 220 K = KI, N
                     EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
     $                      ABS( VL( K, IS+1 ) ) )
  220             CONTINUE
                  REMAX = ONE / EMAX
                  CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
                  CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span>                  DO 230 K = 1, KI - 1
                     VL( K, IS ) = ZERO
                     VL( K, IS+1 ) = ZERO
  230             CONTINUE
               ELSE
                  IF( KI.LT.N-1 ) THEN
                     CALL DGEMV( <span class="string">'N'</span>, N, N-KI-1, ONE, VL( 1, KI+2 ),
     $                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
     $                           VL( 1, KI ), 1 )
                     CALL DGEMV( <span class="string">'N'</span>, N, N-KI-1, ONE, VL( 1, KI+2 ),
     $                           LDVL, WORK( KI+2+N2 ), 1,
     $                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
                  ELSE
                     CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
                     CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  EMAX = ZERO
                  DO 240 K = 1, N
                     EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
     $                      ABS( VL( K, KI+1 ) ) )
  240             CONTINUE
                  REMAX = ONE / EMAX
                  CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
                  CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span>               END IF
<span class="comment">*</span><span class="comment">
</span>            END IF
<span class="comment">*</span><span class="comment">
</span>            IS = IS + 1
            IF( IP.NE.0 )
     $         IS = IS + 1
  250       CONTINUE
            IF( IP.EQ.-1 )
     $         IP = 0
            IF( IP.EQ.1 )
     $         IP = -1
<span class="comment">*</span><span class="comment">
</span>  260    CONTINUE
<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="DTREVC.978"></a><a href="dtrevc.f.html#DTREVC.1">DTREVC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>