AB13ED

Estimating the distance from a real matrix to the nearest complex matrix with an eigenvalue on the imaginary axis, using bisection

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

  To estimate beta(A), the 2-norm distance from a real matrix A to
  the nearest complex matrix with an eigenvalue on the imaginary
  axis. The estimate is given as

         LOW <= beta(A) <= HIGH,

  where either

         (1 + TOL) * LOW >= HIGH,

  or

         LOW = 0   and   HIGH = delta,

  and delta is a small number approximately equal to the square root
  of machine precision times the Frobenius norm (Euclidean norm)
  of A. If A is stable in the sense that all eigenvalues of A lie
  in the open left half complex plane, then beta(A) is the distance
  to the nearest unstable complex matrix, i.e., the complex
  stability radius.

Specification
      SUBROUTINE AB13ED( N, A, LDA, LOW, HIGH, TOL, DWORK, LDWORK,
     $                   INFO )
C     .. Scalar Arguments ..
      DOUBLE PRECISION  HIGH, LOW, TOL
      INTEGER           INFO, LDA, LDWORK, N
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), DWORK(*)

Arguments

Input/Output Parameters

  N       (input) INTEGER
          The order of the matrix A.  N >= 0.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          The leading N-by-N part of this array must contain the
          matrix A.

  LDA     INTEGER
          The leading dimension of array A.  LDA >= MAX(1,N).

  LOW     (output) DOUBLE PRECISION
          A lower bound for beta(A).

  HIGH    (output) DOUBLE PRECISION
          An upper bound for beta(A).

Tolerances
  TOL     DOUBLE PRECISION
          Specifies the accuracy with which LOW and HIGH approximate
          beta(A). If the user sets TOL to be less than SQRT(EPS),
          where EPS is the machine precision (see LAPACK Library
          Routine DLAMCH), then the tolerance is taken to be
          SQRT(EPS).
          The recommended value is TOL = 9, which gives an estimate
          of beta(A) correct to within an order of magnitude.

Workspace
  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) returns the optimal value
          of LDWORK.

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK >= MAX( 1, 3*N*(N+1) ).
          For optimum performance LDWORK should be larger.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  the QR algorithm (LAPACK Library routine DHSEQR)
                fails to converge; this error is very rare.

Method
  Let beta(A) be the 2-norm distance from a real matrix A to the
  nearest complex matrix with an eigenvalue on the imaginary axis.
  It is known that beta(A) = minimum of the smallest singular
  value of (A - jwI), where I is the identity matrix and j**2 = -1,
  and the minimum is taken over all real w.
  The algorithm computes a lower bound LOW and an upper bound HIGH
  for beta(A) by a bisection method in the following way. Given a
  non-negative real number sigma, the Hamiltonian matrix H(sigma)
  is constructed:

                    |   A      -sigma*I |     | A   G  |
        H(sigma) =  |                   | :=  |        | .
                    | sigma*I    -A'    |     | F  -A' |

  It can be shown [1] that H(sigma) has an eigenvalue whose real
  part is zero if and only if sigma >= beta. Any lower and upper
  bounds on beta(A) can be improved by choosing a number between
  them and checking to see if H(sigma) has an eigenvalue with zero
  real part.  This decision is made by computing the eigenvalues of
  H(sigma) using the square reduced algorithm of Van Loan [2].

References
  [1] Byers, R.
      A bisection method for measuring the distance of a stable
      matrix to the unstable matrices.
      SIAM J. Sci. Stat. Comput., Vol. 9, No. 5, pp. 875-880, 1988.

  [2] Van Loan, C.F.
      A symplectic method for approximating all the eigenvalues of a
      Hamiltonian matrix.
      Linear Algebra and its Applications, Vol 61, 233-251, 1984.

Numerical Aspects
  Due to rounding errors the computed values of LOW and HIGH can be
  proven to satisfy

         LOW - p(n) * sqrt(e) * norm(A) <= beta(A)
  and
         beta(A) <= HIGH + p(n) * sqrt(e) * norm(A),

  where p(n) is a modest polynomial of degree 3, e is the machine
  precision and norm(A) is the Frobenius norm of A, see [1].
  The recommended value for TOL is 9 which gives an estimate of
  beta(A) correct to within an order of magnitude.
  AB13ED requires approximately 38*N**3 flops for TOL = 9.

Further Comments
  None
Example

Program Text

*     AB13ED EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX
      PARAMETER        ( NMAX = 20 )
      INTEGER          LDA
      PARAMETER        ( LDA = NMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = 3*NMAX*( NMAX + 1 ) )
*     .. Local Scalars ..
      INTEGER          I, INFO, J, N
      DOUBLE PRECISION HIGH, LOW, TOL
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK)
*     .. External Subroutines ..
      EXTERNAL         AB13ED, UD01MD
*     ..
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
*     Read N, TOL and next A (row wise).
      READ ( NIN, FMT = * ) N, TOL
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99995 ) N
      ELSE
         DO 10 I = 1, N
            READ ( NIN, FMT = * ) ( A(I,J), J = 1, N )
   10    CONTINUE
*
         WRITE ( NOUT, FMT = 99998 ) N, TOL
         CALL UD01MD( N, N, 5, NOUT, A, LDA, 'Matrix A', INFO )
*
         CALL AB13ED( N, A, LDA, LOW, HIGH, TOL, DWORK, LDWORK, INFO )
         IF ( INFO.EQ.0 ) THEN
            WRITE ( NOUT, FMT = 99997 ) LOW, HIGH
         ELSE
            WRITE ( NOUT, FMT = 99996 ) INFO
         END IF
      END IF
      STOP
*
99999 FORMAT (' AB13ED EXAMPLE PROGRAM RESULTS', /1X)
99998 FORMAT (' N =', I4, 2X, 'TOL =', D10.3)
99997 FORMAT (' LOW  =', D18.11, /' HIGH =', D18.11)
99996 FORMAT (' INFO on exit from AB13ED = ', I2)
99995 FORMAT (/' N is out of range.',/' N = ',I5)
      END
Program Data
AB13ED EXAMPLE PROGRAM DATA
5, 9.0D0
1.0D-01  1.0D-00  0.0D-00  0.0D-00  0.0D-00
0.0D-00  1.0D-01  1.0D-00  0.0D-00  0.0D-00
0.0D-00  0.0D-00  1.0D-01  1.0D-00  0.0D-00
0.0D-00  0.0D-00  0.0D-00  1.0D-01  1.0D-00
0.0D-00  0.0D-00  0.0D-00  0.0D-00  1.0D-01
Program Results
 AB13ED EXAMPLE PROGRAM RESULTS

 N =   5  TOL = 0.900D+01
 Matrix A ( 5X 5)

            1              2              3              4              5
  1    0.1000000D+00  0.1000000D+01  0.0000000D+00  0.0000000D+00  0.0000000D+00
  2    0.0000000D+00  0.1000000D+00  0.1000000D+01  0.0000000D+00  0.0000000D+00
  3    0.0000000D+00  0.0000000D+00  0.1000000D+00  0.1000000D+01  0.0000000D+00
  4    0.0000000D+00  0.0000000D+00  0.0000000D+00  0.1000000D+00  0.1000000D+01
  5    0.0000000D+00  0.0000000D+00  0.0000000D+00  0.0000000D+00  0.1000000D+00
 
 LOW  = 0.20929379255D-05
 HIGH = 0.20793050504D-04

Return to index