## AB13FD

### Computing the distance from a real matrix to the nearest complex matrix with an eigenvalue on the imaginary axis, using SVD

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

Purpose

```  To compute beta(A), the 2-norm distance from a real matrix A to
the nearest complex matrix with an eigenvalue on the imaginary
axis. 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 complex
stability radius, i.e., the distance to the nearest unstable
complex matrix. The value of beta(A) is the minimum of the
smallest singular value of (A - jwI), taken over all real w.
The value of w corresponding to the minimum is also computed.

```
Specification
```      SUBROUTINE AB13FD( N, A, LDA, BETA, OMEGA, TOL, DWORK, LDWORK,
\$                   CWORK, LCWORK, INFO )
C     .. Scalar Arguments ..
INTEGER           INFO, LCWORK, LDA, LDWORK, N
DOUBLE PRECISION  BETA, OMEGA, TOL
C     .. Array Arguments ..
DOUBLE PRECISION  A(LDA,*), DWORK(*)
COMPLEX*16        CWORK(*)

```
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).

BETA    (output) DOUBLE PRECISION
The computed value of beta(A), which actually is an upper
bound.

OMEGA   (output) DOUBLE PRECISION
The value of w such that the smallest singular value of
(A - jwI) equals beta(A).

```
Tolerances
```  TOL     DOUBLE PRECISION
Specifies the accuracy with which beta(A) is to be
calculated. (See the Numerical Aspects section below.)
If the user sets TOL to be less than EPS, where EPS is the
machine precision (see LAPACK Library Routine DLAMCH),
then the tolerance is taken to be EPS.

```
Workspace
```  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
If DWORK(1) is not needed, the first 2*N*N entries of
DWORK may overlay CWORK.

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

CWORK   COMPLEX*16 array, dimension (LCWORK)
On exit, if INFO = 0, CWORK(1) returns the optimal value
of LCWORK.
If CWORK(1) is not needed, the first N*N entries of
CWORK may overlay DWORK.

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

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
= 1:  the routine fails to compute beta(A) within the
specified tolerance. Nevertheless, the returned
value is an upper bound on beta(A);
= 2:  either the QR or SVD algorithm (LAPACK Library
routines DHSEQR, DGESVD or ZGESVD) fails to
converge; this error is very rare.

```
Method
```  AB13FD combines the methods of  and  into a provably
reliable, quadratically convergent algorithm. It uses the simple
bisection strategy of  to find an interval which contains
beta(A), and then switches to the modified bisection strategy of
 which converges quadratically to a minimizer. Note that the
efficiency of the strategy degrades if there are several local
minima that are near or equal the global minimum.

```
References
```   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.

 Boyd, S. and Balakrishnan, K.
A regularity result for the singular values of a transfer
matrix and a quadratically convergent algorithm for computing
its L-infinity norm.
Systems and Control Letters, Vol. 15, pp. 1-7, 1990.

```
Numerical Aspects
```  In the presence of rounding errors, the computed function value
BETA  satisfies

beta(A) <= BETA + epsilon,

BETA/(1+TOL) - delta <= MAX(beta(A), SQRT(2*N*EPS)*norm(A)),

where norm(A) is the Frobenius norm of A,

epsilon = p(N) * EPS * norm(A),
and
delta   = p(N) * SQRT(EPS) * norm(A),

and p(N) is a low degree polynomial. It is recommended to choose
TOL greater than SQRT(EPS). Although rounding errors can cause
AB13FD to fail for smaller values of TOL, nevertheless, it usually
succeeds. Regardless of success or failure, the first inequality
holds.

```
```  None
```
Example

Program Text

```*     AB13FD 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 + 2 ) )
INTEGER          LCWORK
PARAMETER        ( LCWORK = NMAX*( NMAX + 3 ) )
*     .. Local Scalars ..
INTEGER          I, INFO, J, N
DOUBLE PRECISION BETA, OMEGA, TOL
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK)
COMPLEX*16       CWORK(LCWORK)
*     .. External Subroutines ..
EXTERNAL         AB13FD, 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, 'A', INFO )
*
CALL AB13FD( N, A, LDA, BETA, OMEGA, TOL, DWORK, LDWORK, CWORK,
\$                LCWORK, INFO )
*
IF ( INFO.NE.0 )
\$      WRITE ( NOUT, FMT = 99996 ) INFO
WRITE ( NOUT, FMT = 99997 ) BETA, OMEGA
END IF
*
99999 FORMAT (' AB13FD EXAMPLE PROGRAM RESULTS', /1X)
99998 FORMAT (' N =', I2, 3X, 'TOL =', D10.3)
99997 FORMAT (' Stability radius :', D18.11, /
*        ' Minimizing omega :', D18.11)
99996 FORMAT (' INFO on exit from AB13FD = ', I2)
99995 FORMAT (/' N is out of range.',/' N = ',I5)
END
```
Program Data
```AB13FD EXAMPLE PROGRAM DATA
4   0.0D-00   0.0D-00
246.500        242.500        202.500       -197.500
-252.500       -248.500       -207.500        202.500
-302.500       -297.500       -248.500        242.500
-307.500       -302.500       -252.500        246.500
```
Program Results
``` AB13FD EXAMPLE PROGRAM RESULTS

N = 4   TOL = 0.000D+00
A ( 4X 4)

1              2              3              4
1    0.2465000D+03  0.2425000D+03  0.2025000D+03 -0.1975000D+03
2   -0.2525000D+03 -0.2485000D+03 -0.2075000D+03  0.2025000D+03
3   -0.3025000D+03 -0.2975000D+03 -0.2485000D+03  0.2425000D+03
4   -0.3075000D+03 -0.3025000D+03 -0.2525000D+03  0.2465000D+03