          DAREX - A Benchmark Collection for DARE
         =========================================

The FORTRAN 77 subroutine DAREX is designed to generate all examples of
discrete-time algebraic Riccati equations (DARE)

                  T            T               T    -1  T       T
(I)  0 = DR(X) = A X A - X - (A X B + S) (R + B X B)  (B X A + S ) + Q

collected in [1]. Here, A, Q, and X are n-by-n matrices, B and S are
n-by-m, and R is m-by-m. The coefficient matrices Q and R are symmetric
and usually, the required solution matrix X is symmetric, too. The
coefficient matrix Q is often given in factored form as 

                T 
(II)     Q  =  C  Q0 C,

where C is p-by-n and Q0 is p-by-p. This factorization often arises
in (but is not limited to) DARE coming from control theory. Also, if R
is nonsingular, the DARE (I) is equivalent to 

                          T             T             -1 
(III)   0  =  DR(X)  =  AA X AA - X - AA X G (I + X G)  X AA + QQ 

where
                     -1 T                     -1 T
        AA =  A - B R  S ,       QQ =  Q - S R  S , 

and
                 -1  T
(IV)    G  =  B R   B .

The required solution X often has some particular or extremal
properties, for instance in control theory one is usually concerned 
with the "stabilizing" solution of (I), i.e., a solution such that the
"feedback gain matrix" 
                                 -1  T       T
(V)	F  =  A  -  B (R + B X B)  (B X A + S )
 
has all its eigenvalues inside the unit circle.

The benchmark collection [1] consists of examples which can be used for
testing purposes in the construction of new numerical methods to solve
DAREs or as a reference set for the comparison of methods. Although
the presented benchmark examples have a control-theoretic background,
they can be considered as examples for general DARE of the form (I).
    
---------------------------------------------------------------------------

IMPLEMENTATION:
===============

DAREX, the required subroutines SP2SY, SY2SP, and the sample program
EXAMPLE were implemented and documented according to the standards
given in [2]. So far only DOUBLE PRECISION versions of these subroutines
are available.  
The subroutines are based upon the DOUBLE PRECISION versions of the BLAS
and LAPACK [3]. Thus, to run the codes, you have to link the BLAS and
LAPACK libraries. If those are not available on your computer, see the
HELP AND BUGS section in this file.  
A call to DAREX looks like

      CALL DAREX(NO, N, M, P, NPAR, DPARAM, A, LDA, B, LDB, C, LDC,
     1           Q, LDQ, R, LDR, S, LDS, X, LDX, NOTE, STORE, WITHC, 
     2           WITHG, WITHS, RWORK, IERR)

A detailed description of all parameters in the argument list of DAREX    
is given in the prolog of the subroutine DAREX.

The coefficient matrices G and Q from (II) and (IV) can be returned either
in the "product form" (i.e., the matrices G and Q are returned in the
arrays R and Q) or in the "factored forms" as in (II) and (IV)
respectively (i.e., the array R contains the matrix R, Q0 is returned in Q
and the matrices B and C are returned in the arrays B, C, respectively) or
in either combination of these forms. The symmetric matrices G, Q, R, Q0
can be stored by columns in a full two-dimensional array (standard matrix
storage mode in FORTRAN 77) or in lower (upper) packed storage mode, i.e.,
the lower (upper) triangle of the matrices is stored by columns.
If an analytical solution of the DARE is available, it will be returned,
too (in standard storage mode) in the array X. This returned solution is
generally the "stabilizing" solution in the control-theoretic sense, i.e.,
the feedback gain matrix F in (V) has all its eigenvalues inside the unit
circle. 

Most of the examples have one INTEGER and/or one or several real
(implemented as DOUBLE PRECISION) parameters. These can be supplied by
the user or they are set by default.  
Some examples require data files which are supplied together with
darex.f. NOTE that the names of the supplied data files (see next
section) contain only capital letters which is obligatory.
A sample program and sample input/output files are also provided (in a
subdirectory named EXAMPLES). Those will be briefly described in the
Section EXAMPLES below.

---------------------------------------------------------------------------

CONTENTS:
=========

You can receive the file darex_f.tar.Z by anonymous ftp at

	ftp.tu-chemnitz.de

from the directory

	pub/Local/mathematik/Benner

(observe the capital "L" in Local !) where you can also receive a
postscript version of the preprint [1] (this file is called blm2.ps.Z).

Then by decompressing darex_f.tar.Z via

	uncompress darex_f.tar.Z

or

	compress -d darex_f.tar.Z

and extracting the resulting file darex.tar by

	tar xf darex_f.tar

a directory darex_f is created containing the following files and
subdirectories: 

EXAMPLES    - A directory containing a sample program, example.f, and a
              sample Makefile as well as ASCII files containing sample
              inputs to and sample outputs from EXAMPLE. Using the sample
              program is described below by showing some sample inputs
              and outputs.
README.f77  - This file.
SOURCE      - A directory containing the source codes required for the
              benchmark collection.  

The subdirectory SOURCE contains the following files:

darex.f     - The FORTRAN 77 subroutine DAREX for generating all the
              benchmark examples presented in [1].
DAREX6.DAT  - Data file in ASCII format required by DAREX for generating
              Example 6 of [1].    
DAREX7.DAT  - Data file in ASCII format required by DAREX for generating
              Example 7 of [1].    
DAREX8.DAT  - Data file in ASCII format required by DAREX for generating
              Example 8 of [1].    
DAREX9.DAT  - Data file in ASCII format required by DAREX for generating
              Example 9 of [1].    
DAREX11.DAT - Data file in ASCII format required by DAREX for generating
              Example 11 of [1].    
sp2sy.f     - The FORTRAN 77 subroutine SP2SY which converts a symmetric 
              matrix stored in packed storage mode to full (standard)
              storage mode. 
sy2sp.f     - The FORTRAN 77 subroutine SY2SP which converts a symmetric 
              matrix stored in full (standard) storage mode to packed
              storage mode. 

---------------------------------------------------------------------------

EXAMPLES:
=========

In the following we will describe how to use DAREX by showing some input
to the sample program EXAMPLE and the resulting output. Note that in
EXAMPLE, R is declared as one-dimensional array designed for packed
storage mode whereas Q is declared as a standard two-dimensional array. 
This is done for demonstration purposes. The actual declaration of Q and 
R should depend upon the desired storage mode (although each declaration
can be used with each storage mode - but this makes life more
complicated).  
In the sample program EXAMPLE, the input is expected from the standard
input device and the output is sent to the standard output device, i.e., 
in both cases the terminal.
After creating an executable file via, e.g., (note that you have to adept
the Makefile to your computer environment)

        make example

there will be an executable file named "darex". You can now use the data 
files provided in the subdirectory EXAMPLES called 

        darex#.d

where # stands for any one of the numbers 1, 4, 14, 15. These numbers
correspond to the example numbers in [1] and if in the following, a #
appears in a filename, this will replace any of those four numbers. (Of
course, you can create data files yourself and give them any name you
want.) 
"darex" can now read the data from such a file by typing

        example < darex#.d

The output files in subdirectory EXAMPLES called

        darex#.r

were created by

        example < darex#.d > darex#.r

If for any reason, "piping" via ">" and "<" is not possible in your
computer environment, you can insert the line

        OPEN(NIN,FILE='darex#.d')

at the beginning of the executable statements in example.f. Then you have
to insert the line 

        CLOSE(NIN)

just before 

        CALL DAREX(....)

or any other place following the last READ directive.
You can create an output file similar by inserting at the beginning of
the executable statements

        OPEN(NOUT,FILE='darex#.r')

and closing the file after the last WRITE command (but before the STOP
directive !) by

        CLOSE(NOUT)

In the following we will assume that "piping" is possible. We now create
some of the benchmark examples by showing some input data (defining the
input paramters to DAREX) for EXAMPLE and the resulting output. 

NOTE that the last line of each input file consists of one character
constant and three boolean expressions. They define the input arguments
STORE, WITHC, WITHG, and WITHS of DAREX. A detailed description of these
parameters is given in the prolog of DAREX in the file darex.f. 

_________________________________________________________________________ 

A) To generate Example 1 of [1] where the arrays R and Q are to be stored
as full matrices (standard 2-dimensional FORTRAN 77 array) (achieved by
STORE = 'F'), Q is returned in factored form, i.e., C is returned and Q
contains the matrix Q0 from (II) (WITHC = .TRUE.), B and R are to be
returned (WITHG = .FALSE.), and S is to be returned (WITHS = .TRUE.),
we need the following input file, darex1.d : 

DAREX EXAMPLE PROGRAM DATA
1 0
'F' .TRUE. .FALSE. .TRUE. 

The first line after the heading contains the number of the required
example (1) and the number of the user-defined parameters (0).
As mentioned above, the next line defines in which form the coefficient
matrices of (I) are returned.
Then by typing

	darex < darex1.d > darex1.r

we obtain the following output file darex1.r :

 DAREX EXAMPLE PROGRAM RESULTS
 
 Laub 1979, Ex. 2: uncontrollable-unobservable data                    

 Order of matrix A:               N  =   2
 Number of columns in matrix B:   M  =   1
 Number of rows in matrix C:      P  =   2

 A  = 
   4.0000   3.0000 
  -4.5000  -3.5000 

 B  = 
   1.0000 
  -1.0000 

 R = 
   1.0000 

 C  = 
   1.0000    .0000 
    .0000   1.0000 

 Q0 = 
   9.0000   6.0000 
   6.0000   4.0000 

 S  = 
    .0000 
    .0000 

 The solution matrix X is 
  14.5623   9.7082 
   9.7082   6.4721 

________________________________________________________________________

B) Now we want to generate Example 4 of [1]. Again we have no input
parameters, but the arrays Q and R are to be stored in lower packed
storage mode. On return we want Q to be formed explicitly (achieved by 
WITHC = .FALSE.), and as before, G is not computed and S is to be
returned. 

The input file darex4.d is then 

DAREX EXAMPLE PROGRAM DATA
4 0
'L' .FALSE. .FALSE. .TRUE. 

From this, we obtain the output file darex4.r :

 DAREX EXAMPLE PROGRAM RESULTS
 
 Ionescu/Weiss 1992: singular R matrix, nonzero S matrix               

 Order of matrix A:               N  =   2
 Number of columns in matrix B:   M  =   2
 Number of rows in matrix C:      P  =   2

 A  = 
    .0000   1.0000 
    .0000  -1.0000 

 B  = 
   1.0000    .0000 
   2.0000   1.0000 

 R = 
   9.0000   3.0000 
   3.0000   1.0000 

 Q  = 
   -.3636   -.3636 
   -.3636    .6364 

 S  = 
   3.0000   1.0000 
  -1.0000   7.0000 

_______________________________________________________________________

C) Next, we create the coefficient matrices corresponding to Example 14
from [1]. In this example there are four real parameters to choose. We
want to set three of them whereas the fourth is to be set by
default. Thus, the first line contains again the number of the example
corresponding to [1], i.e., NO = 14, and the number of user-supplied
parameters, i.e., NPAR = 3. In the next line, the parameters are defined. 
The first one is set to 1000.0, the second is 3.0, and the third is 5.0.

Using the data file darex14.d,

DAREX EXAMPLE PROGRAM DATA
14 3
.1D4 .3D1 .5D1
'U' .TRUE. .FALSE. .FALSE. 
 
we obtain by

	darex < darex14.d > darex14.r

the output file

 DAREX EXAMPLE PROGRAM RESULTS
 
 Pappas et al. 1980: process control of paper machine                  

 Order of matrix A:               N  =   4
 Number of columns in matrix B:   M  =   1
 Number of rows in matrix C:      P  =   1

 A  = 
    .9970    .0000    .0000    .0000 
   1.0000    .0000    .0000    .0000 
    .0000   1.0000    .0000    .0000 
    .0000    .0000   1.0000    .0000 

 B  = 
    .0150 
    .0000 
    .0000 
    .0000 

 R = 
    .2500 

 C  = 
    .0000    .0000    .0000   1.0000 

 Q0 = 
   1.0000 

 The solution matrix X is 
  30.6248    .0000    .0000    .0000 
    .0000   1.0000    .0000    .0000 
    .0000    .0000   1.0000    .0000 
    .0000    .0000    .0000   1.0000 


Here, S is not returned since WITHS = .FALSE. . The fourth parameter in
this example defines the 1-by-1 matrix R. From the output file we see that
on default, R = 0.25.

__________________________________________________________________________

D) The last example given here has the number 15 in [1]. This example has 
two parameters, the order N of the DARE (I), and one DOUBLE PRECISION
parameter r defining the coefficient matrix R which is 1-by-1 as in
Example 14. Here, we want to supply both parameters. N is set to 5 and
R = 8.0. Symmetric matrices are to be stored by their lower triangles and
Q is to be formed explicitly. This time, we want also G to be formed as in
(IV) which is achieved by WITHG = .TRUE. Since S is zero in this example,
we do not want it to be formed and thus set WITHS = .FALSE.

The corresponding data file darex15.d is 

DAREX EXAMPLE PROGRAM DATA
15 2
5
.8D1
'L' .FALSE. .TRUE. .FALSE. 

Here, the user-defined parameters are given in two lines, one
containing the information about the order of the Problem (N = 5), the
other one containing the DOUBLE PRECISION parameter r (= 8.0). 

Generating the output file darex15.r as before, we obtain

 DAREX EXAMPLE PROGRAM RESULTS
 
 Pappas et al. 1980, Ex. 3                                             

 Order of matrix A:               N  =   5
 Number of columns in matrix B:   M  =   1
 Number of rows in matrix C:      P  =   5

 A  = 
    .0000   1.0000    .0000    .0000    .0000 
    .0000    .0000   1.0000    .0000    .0000 
    .0000    .0000    .0000   1.0000    .0000 
    .0000    .0000    .0000    .0000   1.0000 
    .0000    .0000    .0000    .0000    .0000 

 G  = 
    .0000    .0000    .0000    .0000    .0000 
    .0000    .0000    .0000    .0000    .0000 
    .0000    .0000    .0000    .0000    .0000 
    .0000    .0000    .0000    .0000    .0000 
    .0000    .0000    .0000    .0000    .1250 

 Q  = 
   1.0000    .0000    .0000    .0000    .0000 
    .0000   1.0000    .0000    .0000    .0000 
    .0000    .0000   1.0000    .0000    .0000 
    .0000    .0000    .0000   1.0000    .0000 
    .0000    .0000    .0000    .0000   1.0000 

 The solution matrix X is 
   1.0000    .0000    .0000    .0000    .0000 
    .0000   2.0000    .0000    .0000    .0000 
    .0000    .0000   3.0000    .0000    .0000 
    .0000    .0000    .0000   4.0000    .0000 
    .0000    .0000    .0000    .0000   5.0000 

---------------------------------------------------------------------------

HELP AND BUGS:
==============

If you have trouble either compiling or running the codes or if you find
any bugs please send an e-mail message reporting the problem or bug to 

        benner@mathematik.tu-chemnitz.de

We will get in touch with you as soon as possible. We will also appreciate
any proposal for examples to be included in further releases of DAREX.

---------------------------------------------------------------------------

REFERENCES:
===========

[1] P. Benner, A.J. Laub and V.Mehrmann
    'A Collection of Benchmark Examples for the Numerical Solution of
     Algebraic Riccati Equations II: Discrete-Time Case',
    Technical Report SPC 95_23, TU Chemnitz-Zwickau (FRG), 1995.

[2] Numerical Algorithms Group,
    'Implementation and Documentation Standards for the Subroutine 
     Library in Control and Systems Theory SLICOT',
    NAG Publication NP2032, Eindhoven/Oxford, 1990.

[3] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz,
    A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov and D. Sorensen
    'LAPACK Users' Guide', 2nd edition, 
    SIAM, Philadelphia, PA, 1994.

---------------------------------------------------------------------------
 
CONTRIBUTORS:
=============

Peter Benner and Volker Mehrmann
Fakultaet fuer Mathematik
Technische Universitaet Chemnitz-Zwickau (FRG)
e-mail: benner@mathematik.tu-chemnitz.de
        mehrmann@mathematik.tu-chemnitz.de

Alan J. Laub
Department of Electrical and Computer Engineering
University of California
Santa Barbara, CA 93106-9560 (USA)
e-mail: laub@ece.ucsb.edu

---------------------------------------------------------------------------

Peter Benner, February 1, 1995.
