SUBROUTINE SB02RD( JOB, DICO, HINV, TRANA, UPLO, SCAL, SORT, FACT,
     $                   LYAPUN, N, A, LDA, T, LDT, V, LDV, G, LDG, Q,
     $                   LDQ, X, LDX, SEP, RCOND, FERR, WR, WI, S, LDS,
     $                   IWORK, DWORK, LDWORK, BWORK, INFO )
C
C     SLICOT RELEASE 5.0.
C
C     Copyright (c) 2002-2009 NICONET e.V.
C
C     This program is free software: you can redistribute it and/or
C     modify it under the terms of the GNU General Public License as
C     published by the Free Software Foundation, either version 2 of
C     the License, or (at your option) any later version.
C
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with this program.  If not, see
C     <http://www.gnu.org/licenses/>.
C
C     PURPOSE
C
C     To solve for X either the continuous-time algebraic Riccati
C     equation
C                                          -1
C        Q + op(A)'*X + X*op(A) - X*op(B)*R  op(B)'*X = 0,           (1)
C
C     or the discrete-time algebraic Riccati equation
C                                                                -1
C        X = op(A)'*X*op(A) - op(A)'*X*op(B)*(R + op(B)'*X*op(B))  *
C                             op(B)'*X*op(A) + Q,                    (2)
C
C     where op(M) = M or M' (M**T), A, op(B), Q, and R are N-by-N,
C     N-by-M, N-by-N, and M-by-M matrices respectively, with Q symmetric
C     and R symmetric nonsingular; X is an N-by-N symmetric matrix.
C                           -1
C     The matrix G = op(B)*R  *op(B)' must be provided on input, instead
C     of B and R, that is, the continuous-time equation
C
C        Q + op(A)'*X + X*op(A) - X*G*X = 0,                         (3)
C
C     or the discrete-time equation
C                                -1
C        Q + op(A)'*X*(I_n + G*X)  *op(A) - X = 0,                   (4)
C
C     are solved, where G is an N-by-N symmetric matrix. SLICOT Library
C     routine SB02MT should be used to compute G, given B and R. SB02MT
C     also enables to solve Riccati equations corresponding to optimal
C     problems with coupling terms.
C
C     The routine also returns the computed values of the closed-loop
C     spectrum of the optimal system, i.e., the stable eigenvalues
C     lambda(1),...,lambda(N) of the corresponding Hamiltonian or
C     symplectic matrix associated to the optimal problem. It is assumed
C     that the matrices A, G, and Q are such that the associated
C     Hamiltonian or symplectic matrix has N stable eigenvalues, i.e.,
C     with negative real parts, in the continuous-time case, and with
C     moduli less than one, in the discrete-time case.
C
C     Optionally, estimates of the conditioning and error bound on the
C     solution of the Riccati equation (3) or (4) are returned.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOB     CHARACTER*1
C             Specifies the computation to be performed, as follows:
C             = 'X':  Compute the solution only;
C             = 'C':  Compute the reciprocal condition number only;
C             = 'E':  Compute the error bound only;
C             = 'A':  Compute all: the solution, reciprocal condition
C                     number, and the error bound.
C
C     DICO    CHARACTER*1
C             Specifies the type of Riccati equation to be solved or
C             analyzed, as follows:
C             = 'C':  Equation (3), continuous-time case;
C             = 'D':  Equation (4), discrete-time case.
C
C     HINV    CHARACTER*1
C             If DICO = 'D' and JOB = 'X' or JOB = 'A', specifies which
C             symplectic matrix is to be constructed, as follows:
C             = 'D':  The matrix H in (6) (see METHOD) is constructed;
C             = 'I':  The inverse of the matrix H in (6) is constructed.
C             HINV is not used if DICO = 'C', or JOB = 'C' or 'E'.
C
C     TRANA   CHARACTER*1
C             Specifies the form of op(A) to be used, as follows:
C             = 'N':  op(A) = A    (No transpose);
C             = 'T':  op(A) = A**T (Transpose);
C             = 'C':  op(A) = A**T (Conjugate transpose = Transpose).
C
C     UPLO    CHARACTER*1
C             Specifies which triangle of the matrices G and Q is
C             stored, as follows:
C             = 'U':  Upper triangle is stored;
C             = 'L':  Lower triangle is stored.
C
C     SCAL    CHARACTER*1
C             If JOB = 'X' or JOB = 'A', specifies whether or not a
C             scaling strategy should be used, as follows:
C             = 'G':  General scaling should be used;
C             = 'N':  No scaling should be used.
C             SCAL is not used if JOB = 'C' or 'E'.
C
C     SORT    CHARACTER*1
C             If JOB = 'X' or JOB = 'A', specifies which eigenvalues
C             should be obtained in the top of the Schur form, as
C             follows:
C             = 'S':  Stable   eigenvalues come first;
C             = 'U':  Unstable eigenvalues come first.
C             SORT is not used if JOB = 'C' or 'E'.
C
C     FACT    CHARACTER*1
C             If JOB <> 'X', specifies whether or not a real Schur
C             factorization of the closed-loop system matrix Ac is
C             supplied on entry, as follows:
C             = 'F':  On entry, T and V contain the factors from a real
C                     Schur factorization of the matrix Ac;
C             = 'N':  A Schur factorization of Ac will be computed
C                     and the factors will be stored in T and V.
C             For a continuous-time system, the matrix Ac is given by
C                Ac = A - G*X, if TRANA = 'N', or
C                Ac = A - X*G, if TRANA = 'T' or 'C',
C             and for a discrete-time system, the matrix Ac is given by
C                Ac = inv(I_n + G*X)*A, if TRANA = 'N', or
C                Ac = A*inv(I_n + X*G), if TRANA = 'T' or 'C'.
C             FACT is not used if JOB = 'X'.
C
C     LYAPUN  CHARACTER*1
C             If JOB <> 'X', specifies whether or not the original or
C             "reduced" Lyapunov equations should be solved for
C             estimating reciprocal condition number and/or the error
C             bound, as follows:
C             = 'O':  Solve the original Lyapunov equations, updating
C                     the right-hand sides and solutions with the
C                     matrix V, e.g., X <-- V'*X*V;
C             = 'R':  Solve reduced Lyapunov equations only, without
C                     updating the right-hand sides and solutions.
C                     This means that a real Schur form T of Ac appears
C                     in the equations, instead of Ac.
C             LYAPUN is not used if JOB = 'X'.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrices A, Q, G, and X.  N >= 0.
C
C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
C             If JOB = 'X' or JOB = 'A' or FACT = 'N' or LYAPUN = 'O',
C             the leading N-by-N part of this array must contain the
C             coefficient matrix A of the equation.
C             If JOB = 'C' or 'E' and FACT = 'F' and LYAPUN = 'R', A is
C             not referenced.
C
C     LDA     INTEGER
C             The leading dimension of the array A.
C             LDA >= MAX(1,N), if JOB  = 'X' or JOB = 'A' or
C                                 FACT = 'N' or LYAPUN = 'O'.
C             LDA >= 1,        otherwise.
C
C     T       (input or output) DOUBLE PRECISION array, dimension
C             (LDT,N)
C             If JOB <> 'X' and FACT = 'F', then T is an input argument
C             and on entry, the leading N-by-N upper Hessenberg part of
C             this array must contain the upper quasi-triangular matrix
C             T in Schur canonical form from a Schur factorization of Ac
C             (see argument FACT).
C             If JOB <> 'X' and FACT = 'N', then T is an output argument
C             and on exit, if INFO = 0 or INFO = 7, the leading N-by-N
C             upper Hessenberg part of this array contains the upper
C             quasi-triangular matrix T in Schur canonical form from a
C             Schur factorization of Ac (see argument FACT).
C             If JOB = 'X', the array T is not referenced.
C
C     LDT     INTEGER
C             The leading dimension of the array T.
C             LDT >= 1,        if JOB =  'X';
C             LDT >= MAX(1,N), if JOB <> 'X'.
C
C     V       (input or output) DOUBLE PRECISION array, dimension
C             (LDV,N)
C             If JOB <> 'X' and FACT = 'F', then V is an input argument
C             and on entry, the leading N-by-N part of this array must
C             contain the orthogonal matrix V from a real Schur
C             factorization of Ac (see argument FACT).
C             If JOB <> 'X' and FACT = 'N', then V is an output argument
C             and on exit, if INFO = 0 or INFO = 7, the leading N-by-N
C             part of this array contains the orthogonal N-by-N matrix
C             from a real Schur factorization of Ac (see argument FACT).
C             If JOB = 'X', the array V is not referenced.
C
C     LDV     INTEGER
C             The leading dimension of the array V.
C             LDV >= 1,        if JOB =  'X';
C             LDV >= MAX(1,N), if JOB <> 'X'.
C
C     G       (input/output) DOUBLE PRECISION array, dimension (LDG,N)
C             On entry, the leading N-by-N upper triangular part (if
C             UPLO = 'U') or lower triangular part (if UPLO = 'L') of
C             this array must contain the upper triangular part or lower
C             triangular part, respectively, of the symmetric matrix G.
C             On exit, if JOB = 'X' and DICO = 'D', or JOB <> 'X' and
C             LYAPUN = 'R', the leading N-by-N part of this array
C             contains the symmetric matrix G fully stored.
C             If JOB <> 'X' and LYAPUN = 'R', this array is modified
C             internally, but restored on exit.
C
C     LDG     INTEGER
C             The leading dimension of the array G.  LDG >= MAX(1,N).
C
C     Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
C             On entry, the leading N-by-N upper triangular part (if
C             UPLO = 'U') or lower triangular part (if UPLO = 'L') of
C             this array must contain the upper triangular part or lower
C             triangular part, respectively, of the symmetric matrix Q.
C             On exit, if JOB = 'X' and DICO = 'D', or JOB <> 'X' and
C             LYAPUN = 'R', the leading N-by-N part of this array
C             contains the symmetric matrix Q fully stored.
C             If JOB <> 'X' and LYAPUN = 'R', this array is modified
C             internally, but restored on exit.
C
C     LDQ     INTEGER
C             The leading dimension of the array Q.  LDQ >= MAX(1,N).
C
C     X       (input or output) DOUBLE PRECISION array, dimension
C             (LDX,N)
C             If JOB = 'C' or JOB = 'E', then X is an input argument
C             and on entry, the leading N-by-N part of this array must
C             contain the symmetric solution matrix of the algebraic
C             Riccati equation. If LYAPUN = 'R', this array is modified
C             internally, but restored on exit; however, it could differ
C             from the input matrix at the round-off error level.
C             If JOB = 'X' or JOB = 'A', then X is an output argument
C             and on exit, if INFO = 0 or INFO >= 6, the leading N-by-N
C             part of this array contains the symmetric solution matrix
C             X of the algebraic Riccati equation.
C
C     LDX     INTEGER
C             The leading dimension of the array X.  LDX >= MAX(1,N).
C
C     SEP     (output) DOUBLE PRECISION
C             If JOB = 'C' or JOB = 'A', and INFO = 0 or INFO = 7, the
C             estimated quantity
C                sep(op(Ac),-op(Ac)'), if DICO = 'C', or
C                sepd(op(Ac),op(Ac)'), if DICO = 'D'. (See METHOD.)
C             If JOB = 'C' or JOB = 'A' and X = 0, or JOB = 'E', SEP is
C             not referenced.
C             If JOB = 'X', and INFO = 0, INFO = 5 or INFO = 7,
C             SEP contains the scaling factor used, which should
C             multiply the (2,1) submatrix of U to recover X from the
C             first N columns of U (see METHOD). If SCAL = 'N', SEP is
C             set to 1.
C
C     RCOND   (output) DOUBLE PRECISION
C             If JOB = 'C' or JOB = 'A', and INFO = 0 or INFO = 7, an
C             estimate of the reciprocal condition number of the
C             algebraic Riccati equation.
C             If N = 0 or X = 0, RCOND is set to 1 or 0, respectively.
C             If JOB = 'X', or JOB = 'E', RCOND is not referenced.
C
C     FERR    (output) DOUBLE PRECISION
C             If JOB = 'E' or JOB = 'A', and INFO = 0 or INFO = 7, an
C             estimated forward error bound for the solution X. If XTRUE
C             is the true solution, FERR bounds the magnitude of the
C             largest entry in (X - XTRUE) divided by the magnitude of
C             the largest entry in X.
C             If N = 0 or X = 0, FERR is set to 0.
C             If JOB = 'X', or JOB = 'C', FERR is not referenced.
C
C     WR      (output) DOUBLE PRECISION array, dimension (2*N)
C     WI      (output) DOUBLE PRECISION array, dimension (2*N)
C             If JOB = 'X' or JOB = 'A', and INFO = 0 or INFO >= 5,
C             these arrays contain the real and imaginary parts,
C             respectively, of the eigenvalues of the 2N-by-2N matrix S,
C             ordered as specified by SORT (except for the case
C             HINV = 'D', when the order is opposite to that specified
C             by SORT). The leading N elements of these arrays contain
C             the closed-loop spectrum of the system matrix Ac (see
C             argument FACT). Specifically,
C                lambda(k) = WR(k) + j*WI(k), for k = 1,2,...,N.
C             If JOB = 'C' or JOB = 'E', these arrays are not
C             referenced.
C
C     S       (output) DOUBLE PRECISION array, dimension (LDS,2*N)
C             If JOB = 'X' or JOB = 'A', and INFO = 0 or INFO >= 5, the
C             leading 2N-by-2N part of this array contains the ordered
C             real Schur form S of the (scaled, if SCAL = 'G')
C             Hamiltonian or symplectic matrix H. That is,
C
C                    ( S    S   )
C                    (  11   12 )
C                S = (          ),
C                    ( 0    S   )
C                    (       22 )
C
C             where S  , S   and S   are N-by-N matrices.
C                    11   12      22
C             If JOB = 'C' or JOB = 'E', this array is not referenced.
C
C     LDS     INTEGER
C             The leading dimension of the array S.
C             LDS >= MAX(1,2*N), if JOB = 'X' or JOB = 'A';
C             LDS >= 1,          if JOB = 'C' or JOB = 'E'.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (LIWORK)
C             LIWORK >= 2*N,          if JOB = 'X';
C             LIWORK >= N*N,          if JOB = 'C' or JOB = 'E';
C             LIWORK >= MAX(2*N,N*N), if JOB = 'A'.
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, or INFO = 7, DWORK(1) returns the
C             optimal value of LDWORK. If INFO = 0, or INFO >= 5, and
C             JOB = 'X', or JOB = 'A', then DWORK(2) returns an estimate
C             RCONDU of the reciprocal of the condition number (in the
C             1-norm) of the N-th order system of algebraic equations
C             from which the solution matrix X is obtained, and DWORK(3)
C             returns the reciprocal pivot growth factor for the LU
C             factorization of the coefficient matrix of that system
C             (see SLICOT Library routine MB02PD); if DWORK(3) is much
C             less than 1, then the computed X and RCONDU could be
C             unreliable.
C             If DICO = 'D', and JOB = 'X', or JOB = 'A', then DWORK(4)
C             returns the reciprocal condition number RCONDA of the
C             given matrix A, and DWORK(5) returns the reciprocal pivot
C             growth factor for A or for its leading columns, if A is
C             singular (see SLICOT Library routine MB02PD); if DWORK(5)
C             is much less than 1, then the computed S and RCONDA could
C             be unreliable.
C             On exit, if INFO = 0, or INFO >= 4, and JOB = 'X', the
C             elements DWORK(6:5+4*N*N) contain the 2*N-by-2*N
C             transformation matrix  U  which reduced the Hamiltonian or
C             symplectic matrix  H  to the ordered real Schur form  S.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= 5+MAX(1,4*N*N+8*N), if JOB = 'X' or JOB = 'A';
C             This may also be used for JOB = 'C' or JOB = 'E', but
C             exact bounds are as follows:
C             LDWORK >= 5 + MAX(1,LWS,LWE) + LWN, where
C             LWS = 0,       if FACT = 'F' or  LYAPUN = 'R';
C                 = 5*N,     if FACT = 'N' and LYAPUN = 'O' and
C                                              DICO = 'C' and JOB = 'C';
C                 = 5*N+N*N, if FACT = 'N' and LYAPUN = 'O' and
C                                              DICO = 'C' and JOB = 'E';
C                 = 5*N+N*N, if FACT = 'N' and LYAPUN = 'O' and
C                                              DICO = 'D';
C             LWE = 2*N*N,                if DICO = 'C' and JOB = 'C';
C                 = 4*N*N,                if DICO = 'C' and JOB = 'E';
C                 = MAX(3,2*N*N) + N*N,   if DICO = 'D' and JOB = 'C';
C                 = MAX(3,2*N*N) + 2*N*N, if DICO = 'D' and JOB = 'E';
C             LWN = 0,   if LYAPUN = 'O' or   JOB = 'C';
C                 = 2*N, if LYAPUN = 'R' and DICO = 'C' and JOB = 'E';
C                 = 3*N, if LYAPUN = 'R' and DICO = 'D' and JOB = 'E'.
C             For optimum performance LDWORK should sometimes be larger.
C
C     BWORK   LOGICAL array, dimension (LBWORK)
C             LBWORK >= 2*N,          if JOB = 'X' or JOB = 'A';
C             LBWORK >= 1,            if JOB = 'C' or JOB = 'E', and
C                                     FACT = 'N' and LYAPUN = 'R';
C             LBWORK >= 0,            otherwise.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             = 1:  if matrix A is (numerically) singular in discrete-
C                   time case;
C             = 2:  if the Hamiltonian or symplectic matrix H cannot be
C                   reduced to real Schur form;
C             = 3:  if the real Schur form of the Hamiltonian or
C                   symplectic matrix H cannot be appropriately ordered;
C             = 4:  if the Hamiltonian or symplectic matrix H has less
C                   than N stable eigenvalues;
C             = 5:  if the N-th order system of linear algebraic
C                   equations, from which the solution matrix X would
C                   be obtained, is singular to working precision;
C             = 6:  if the QR algorithm failed to complete the reduction
C                   of the matrix Ac to Schur canonical form, T;
C             = 7:  if T and -T' have some almost equal eigenvalues, if
C                   DICO = 'C', or T has almost reciprocal eigenvalues,
C                   if DICO = 'D'; perturbed values were used to solve
C                   Lyapunov equations, but the matrix T, if given (for
C                   FACT = 'F'), is unchanged. (This is a warning
C                   indicator.)
C
C     METHOD
C
C     The method used is the Schur vector approach proposed by Laub [1],
C     but with an optional scaling, which enhances the numerical
C     stability [6]. It is assumed that [A,B] is a stabilizable pair
C     (where for (3) or (4), B is any matrix such that B*B' = G with
C     rank(B) = rank(G)), and [E,A] is a detectable pair, where E is any
C     matrix such that E*E' = Q with rank(E) = rank(Q). Under these
C     assumptions, any of the algebraic Riccati equations (1)-(4) is
C     known to have a unique non-negative definite solution. See [2].
C     Now consider the 2N-by-2N Hamiltonian or symplectic matrix
C
C                 ( op(A)   -G    )
C            H =  (               ),                                 (5)
C                 (  -Q   -op(A)' ),
C
C     for continuous-time equation, and
C                         -1              -1
C                 (  op(A)           op(A)  *G       )
C            H =  (        -1                   -1   ),              (6)
C                 ( Q*op(A)     op(A)' + Q*op(A)  *G )
C
C     for discrete-time equation, respectively, where
C                       -1
C            G = op(B)*R  *op(B)'.
C     The assumptions guarantee that H in (5) has no pure imaginary
C     eigenvalues, and H in (6) has no eigenvalues on the unit circle.
C     If Y is an N-by-N matrix then there exists an orthogonal matrix U
C     such that U'*Y*U is an upper quasi-triangular matrix. Moreover, U
C     can be chosen so that the 2-by-2 and 1-by-1 diagonal blocks
C     (corresponding to the complex conjugate eigenvalues and real
C     eigenvalues respectively) appear in any desired order. This is the
C     ordered real Schur form. Thus, we can find an orthogonal
C     similarity transformation U which puts (5) or (6) in ordered real
C     Schur form
C
C            U'*H*U = S = (S(1,1)  S(1,2))
C                         (  0     S(2,2))
C
C     where S(i,j) is an N-by-N matrix and the eigenvalues of S(1,1)
C     have negative real parts in case of (5), or moduli greater than
C     one in case of (6). If U is conformably partitioned into four
C     N-by-N blocks
C
C               U = (U(1,1)  U(1,2))
C                   (U(2,1)  U(2,2))
C
C     with respect to the assumptions we then have
C     (a) U(1,1) is invertible and X = U(2,1)*inv(U(1,1)) solves (1),
C         (2), (3), or (4) with X = X' and non-negative definite;
C     (b) the eigenvalues of S(1,1) (if DICO = 'C') or S(2,2) (if
C         DICO = 'D') are equal to the eigenvalues of optimal system
C         (the 'closed-loop' spectrum).
C
C     [A,B] is stabilizable if there exists a matrix F such that (A-BF)
C     is stable. [E,A] is detectable if [A',E'] is stabilizable.
C
C     The condition number of a Riccati equation is estimated as
C
C     cond = ( norm(Theta)*norm(A) + norm(inv(Omega))*norm(Q) +
C                 norm(Pi)*norm(G) ) / norm(X),
C
C     where Omega, Theta and Pi are linear operators defined by
C
C     Omega(W) = op(Ac)'*W + W*op(Ac),
C     Theta(W) = inv(Omega(op(W)'*X + X*op(W))),
C        Pi(W) = inv(Omega(X*W*X)),
C
C     in the continuous-time case, and
C
C     Omega(W) = op(Ac)'*W*op(Ac) - W,
C     Theta(W) = inv(Omega(op(W)'*X*op(Ac) + op(Ac)'X*op(W))),
C        Pi(W) = inv(Omega(op(Ac)'*X*W*X*op(Ac))),
C
C     in the discrete-time case, and Ac has been defined (see argument
C     FACT). Details are given in the comments of SLICOT Library
C     routines SB02QD and SB02SD.
C
C     The routine estimates the quantities
C
C     sep(op(Ac),-op(Ac)') = 1 / norm(inv(Omega)),
C     sepd(op(Ac),op(Ac)') = 1 / norm(inv(Omega)),
C
C     norm(Theta) and norm(Pi) using 1-norm condition estimator.
C
C     The forward error bound is estimated using a practical error bound
C     similar to the one proposed in [5].
C
C     REFERENCES
C
C     [1] Laub, A.J.
C         A Schur Method for Solving Algebraic Riccati equations.
C         IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
C
C     [2] Wonham, W.M.
C         On a matrix Riccati equation of stochastic control.
C         SIAM J. Contr., 6, pp. 681-697, 1968.
C
C     [3] Sima, V.
C         Algorithms for Linear-Quadratic Optimization.
C         Pure and Applied Mathematics: A Series of Monographs and
C         Textbooks, vol. 200, Marcel Dekker, Inc., New York, 1996.
C
C     [4] Ghavimi, A.R. and Laub, A.J.
C         Backward error, sensitivity, and refinement of computed
C         solutions of algebraic Riccati equations.
C         Numerical Linear Algebra with Applications, vol. 2, pp. 29-49,
C         1995.
C
C     [5] Higham, N.J.
C         Perturbation theory and backward error for AX-XB=C.
C         BIT, vol. 33, pp. 124-136, 1993.
C
C     [6] Petkov, P.Hr., Konstantinov, M.M., and Mehrmann, V.
C         DGRSVX and DMSRIC: Fortran 77 subroutines for solving
C         continuous-time matrix algebraic Riccati equations with
C         condition and accuracy estimates.
C         Preprint SFB393/98-16, Fak. f. Mathematik, Tech. Univ.
C         Chemnitz, May 1998.
C
C     NUMERICAL ASPECTS
C                               3
C     The algorithm requires 0(N ) operations. The solution accuracy
C     can be controlled by the output parameter FERR.
C
C     FURTHER COMMENTS
C
C     To obtain a stabilizing solution of the algebraic Riccati
C     equation for DICO = 'D', set SORT = 'U', if HINV = 'D', or set
C     SORT = 'S', if HINV = 'I'.
C
C     The routine can also compute the anti-stabilizing solutions of
C     the algebraic Riccati equations, by specifying
C         SORT = 'U' if DICO = 'D' and HINV = 'I', or DICO = 'C', or
C         SORT = 'S' if DICO = 'D' and HINV = 'D'.
C
C     Usually, the combinations HINV = 'D' and SORT = 'U', or HINV = 'I'
C     and SORT = 'U', for stabilizing and anti-stabilizing solutions,
C     respectively, will be faster then the other combinations [3].
C
C     The option LYAPUN = 'R' may produce slightly worse or better
C     estimates, and it is faster than the option 'O'.
C
C     This routine is a functionally extended and more accurate
C     version of the SLICOT Library routine SB02MD. Transposed problems
C     can be dealt with as well. Iterative refinement is used whenever
C     useful to solve linear algebraic systems. Condition numbers and
C     error bounds on the solutions are optionally provided.
C
C     CONTRIBUTOR
C
C     V. Sima, Research Institute for Informatics, Bucharest, Apr. 1999.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Oct. 2001,
C     Dec. 2002, Oct. 2004.
C
C     KEYWORDS
C
C     Algebraic Riccati equation, closed loop system, continuous-time
C     system, discrete-time system, optimal regulator, Schur form.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, HALF, ONE
      PARAMETER         ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      CHARACTER         DICO, FACT, HINV, JOB, LYAPUN, SCAL, SORT,
     $                  TRANA, UPLO
      INTEGER           INFO, LDA, LDG, LDQ, LDS, LDT, LDV, LDWORK, LDX,
     $                  N
      DOUBLE PRECISION  FERR, RCOND, SEP
C     .. Array Arguments ..
      LOGICAL           BWORK(*)
      INTEGER           IWORK(*)
      DOUBLE PRECISION  A(LDA,*), DWORK(*), G(LDG,*), Q(LDQ,*),
     $                  S(LDS,*), T(LDT,*), V(LDV,*), WI(*), WR(*),
     $                  X(LDX,*)
C     .. Local Scalars ..
      LOGICAL           COLEQU, DISCR, JBXA, JOBA, JOBC, JOBE, JOBX,
     $                  LHINV, LSCAL, LSCL, LSORT, LUPLO, NOFACT,
     $                  NOTRNA, ROWEQU, UPDATE
      CHARACTER         EQUED, JOBS, LOFACT, LOUP, TRANAT
      INTEGER           I, IERR, IU, IW, IWB, IWC, IWF, IWI, IWR, LDW,
     $                  LWE, LWN, LWS, N2, NN, NP1, NROT
      DOUBLE PRECISION  GNORM, QNORM, PIVOTA, PIVOTU, RCONDA, RCONDU,
     $                  WRKOPT
C     .. External Functions ..
      LOGICAL           LSAME, SB02MR, SB02MS, SB02MV, SB02MW
      DOUBLE PRECISION  DLAMCH, DLANGE, DLANSY
      EXTERNAL          DLAMCH, DLANGE, DLANSY, LSAME, SB02MR, SB02MS,
     $                  SB02MV, SB02MW
C     .. External Subroutines ..
      EXTERNAL          DAXPY, DCOPY, DGEES, DGESV, DLACPY, DLASCL,
     $                  DLASET, DSCAL, DSWAP, DSYMM, MA02AD, MA02ED,
     $                  MB01RU, MB01SD, MB02PD, SB02QD, SB02RU, SB02SD,
     $                  XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         DBLE, MAX
C     .. Executable Statements ..
C
C     Decode the input parameters.
C
      N2  = N + N
      NN  = N*N
      NP1 = N + 1
      INFO = 0
      JOBA   = LSAME( JOB,    'A' )
      JOBC   = LSAME( JOB,    'C' )
      JOBE   = LSAME( JOB,    'E' )
      JOBX   = LSAME( JOB,    'X' )
      NOFACT = LSAME( FACT,   'N' )
      NOTRNA = LSAME( TRANA,  'N' )
      DISCR  = LSAME( DICO,   'D' )
      LUPLO  = LSAME( UPLO,   'U' )
      LSCAL  = LSAME( SCAL,   'G' )
      LSORT  = LSAME( SORT,   'S' )
      UPDATE = LSAME( LYAPUN, 'O' )
      JBXA   = JOBX .OR. JOBA
      LHINV  = .FALSE.
      IF ( DISCR .AND. JBXA )
     $   LHINV = LSAME( HINV, 'D' )
C
C     Test the input scalar arguments.
C
      IF( .NOT.( JBXA .OR. JOBC .OR. JOBE ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( DISCR .OR. LSAME( DICO, 'C' ) ) ) THEN
         INFO = -2
      ELSE IF( DISCR .AND. JBXA ) THEN
         IF( .NOT.( LHINV .OR. LSAME( HINV, 'I' ) ) )
     $      INFO = -3
      END IF
      IF( INFO.EQ.0 ) THEN
         IF( .NOT.( NOTRNA .OR. LSAME( TRANA, 'T' ) .OR.
     $                          LSAME( TRANA, 'C' ) ) ) THEN
            INFO = -4
         ELSE IF( .NOT.( LUPLO .OR. LSAME( UPLO, 'L' ) ) )
     $      THEN
            INFO = -5
         ELSE IF( JBXA ) THEN
            IF( .NOT.( LSCAL .OR. LSAME( SCAL, 'N' ) ) ) THEN
               INFO = -6
            ELSE IF( .NOT.( LSORT .OR. LSAME( SORT, 'U' ) ) ) THEN
               INFO = -7
            END IF
         END IF
         IF( INFO.EQ.0 .AND. .NOT.JOBX ) THEN
            IF( .NOT.( NOFACT .OR. LSAME( FACT, 'F' ) ) ) THEN
               INFO = -8
            ELSE IF( .NOT.( UPDATE .OR. LSAME( LYAPUN, 'R' ) ) ) THEN
               INFO = -9
            END IF
         END IF
         IF( INFO.EQ.0 ) THEN
            IF( N.LT.0 ) THEN
               INFO = -10
            ELSE IF( LDA.LT.1 .OR. ( ( JBXA .OR. NOFACT .OR. UPDATE )
     $         .AND. LDA.LT.N ) ) THEN
               INFO = -12
            ELSE IF( LDT.LT.1 .OR. ( .NOT. JOBX .AND. LDT.LT.N ) ) THEN
               INFO = -14
            ELSE IF( LDV.LT.1 .OR. ( .NOT. JOBX .AND. LDV.LT.N ) ) THEN
               INFO = -16
            ELSE IF( LDG.LT.MAX( 1, N ) ) THEN
               INFO = -18
            ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
               INFO = -20
            ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
               INFO = -22
            ELSE IF( LDS.LT.1 .OR. ( JBXA .AND. LDS.LT.N2 ) ) THEN
               INFO = -29
            ELSE
               IF( JBXA ) THEN
                  IF( LDWORK.LT.5 + MAX( 1, 4*NN + 8*N ) )
     $               INFO = -32
               ELSE
                  IF( NOFACT .AND. UPDATE ) THEN
                     IF( .NOT.DISCR .AND. JOBC ) THEN
                        LWS = 5*N
                     ELSE
                        LWS = 5*N + NN
                     END IF
                  ELSE
                     LWS = 0
                  END IF
                  IF( DISCR ) THEN
                     IF( JOBC ) THEN
                        LWE = MAX( 3, 2*NN) + NN
                     ELSE
                        LWE = MAX( 3, 2*NN) + 2*NN
                     END IF
                  ELSE
                     IF( JOBC ) THEN
                        LWE = 2*NN
                     ELSE
                        LWE = 4*NN
                     END IF
                  END IF
                  IF( UPDATE .OR. JOBC ) THEN
                     LWN = 0
                  ELSE
                     IF( DISCR ) THEN
                        LWN = 3*N
                     ELSE
                        LWN = 2*N
                     END IF
                  END IF
                  IF( LDWORK.LT.5 + MAX( 1, LWS, LWE ) + LWN )
     $               INFO = -32
               END IF
            END IF
         END IF
      END IF
C
      IF ( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'SB02RD', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( N.EQ.0 ) THEN
         IF( JOBX )
     $      SEP = ONE
         IF( JOBC .OR. JOBA )
     $      RCOND = ONE
         IF( JOBE .OR. JOBA )
     $      FERR  = ZERO
         DWORK(1) = ONE
         DWORK(2) = ONE
         DWORK(3) = ONE
         IF ( DISCR ) THEN
            DWORK(4) = ONE
            DWORK(5) = ONE
         END IF
         RETURN
      END IF
C
      IF ( JBXA ) THEN
C
C        Compute the solution matrix X.
C
C        Initialise the Hamiltonian or symplectic matrix associated with
C        the problem.
C        Workspace:  need   0    if DICO = 'C';
C                           6*N, if DICO = 'D'.
C
         CALL SB02RU( DICO, HINV, TRANA, UPLO, N, A, LDA, G, LDG, Q,
     $                LDQ, S, LDS, IWORK, DWORK, LDWORK, IERR )
C
         IF ( IERR.NE.0 ) THEN
            INFO = 1
            IF ( DISCR ) THEN
               DWORK(4) = DWORK(1)
               DWORK(5) = DWORK(2)
            END IF
            RETURN
         END IF
C
         IF ( DISCR ) THEN
            WRKOPT = 6*N
            RCONDA = DWORK(1)
            PIVOTA = DWORK(2)
         ELSE
            WRKOPT = 0
         END IF
C
         IF ( LSCAL ) THEN
C
C           Scale the Hamiltonian or symplectic matrix S, using the
C           square roots of the norms of the matrices Q and G.
C
            QNORM = SQRT( DLANSY( '1-norm', UPLO, N, Q, LDQ, DWORK ) )
            GNORM = SQRT( DLANSY( '1-norm', UPLO, N, G, LDG, DWORK ) )
C
            LSCL = QNORM.GT.GNORM .AND. GNORM.GT.ZERO
            IF( LSCL ) THEN
               CALL DLASCL( 'G', 0, 0, QNORM, GNORM, N, N, S(NP1,1),
     $                      LDS, IERR )
               CALL DLASCL( 'G', 0, 0, GNORM, QNORM, N, N, S(1,NP1),
     $                      LDS, IERR )
            END IF
         ELSE
            LSCL = .FALSE.
         END IF
C
C        Find the ordered Schur factorization of S,  S = U*H*U'.
C        Workspace:  need   5 + 4*N*N + 6*N;
C                    prefer larger.
C
         IU  = 6
         IW  = IU + 4*NN
         LDW = LDWORK - IW + 1
         IF ( .NOT.DISCR ) THEN
            IF ( LSORT ) THEN
               CALL DGEES( 'Vectors', 'Sorted', SB02MV, N2, S, LDS,
     $                     NROT, WR, WI, DWORK(IU), N2, DWORK(IW), LDW,
     $                     BWORK, IERR )
            ELSE
               CALL DGEES( 'Vectors', 'Sorted', SB02MR, N2, S, LDS,
     $                     NROT, WR, WI, DWORK(IU), N2, DWORK(IW), LDW,
     $                     BWORK, IERR )
            END IF
         ELSE
            IF ( LSORT ) THEN
               CALL DGEES( 'Vectors', 'Sorted', SB02MW, N2, S, LDS,
     $                     NROT, WR, WI, DWORK(IU), N2, DWORK(IW), LDW,
     $                     BWORK, IERR )
            ELSE
               CALL DGEES( 'Vectors', 'Sorted', SB02MS, N2, S, LDS,
     $                     NROT, WR, WI, DWORK(IU), N2, DWORK(IW), LDW,
     $                     BWORK, IERR )
            END IF
            IF ( LHINV ) THEN
               CALL DSWAP( N, WR, 1, WR(NP1), 1 )
               CALL DSWAP( N, WI, 1, WI(NP1), 1 )
            END IF
         END IF
         IF ( IERR.GT.N2 ) THEN
            INFO = 3
         ELSE IF ( IERR.GT.0 ) THEN
            INFO = 2
         ELSE IF ( NROT.NE.N ) THEN
            INFO = 4
         END IF
         IF ( INFO.NE.0 ) THEN
            IF ( DISCR ) THEN
               DWORK(4) = RCONDA
               DWORK(5) = PIVOTA
            END IF
            RETURN
         END IF
C
         WRKOPT = MAX( WRKOPT, DWORK(IW) + DBLE( IW - 1 ) )
C
C        Compute the solution of X*U(1,1) = U(2,1) using
C        LU factorization and iterative refinement. The (2,1) block of S
C        is used as a workspace for factoring U(1,1).
C        Workspace:  need   5 + 4*N*N + 8*N.
C
C        First transpose U(2,1) in-situ.
C
         DO 20 I = 1, N - 1
            CALL DSWAP( N-I, DWORK(IU+N+I*(N2+1)-1), N2,
     $                  DWORK(IU+N+(I-1)*(N2+1)+1), 1 )
   20    CONTINUE
C
         IWR = IW
         IWC = IWR + N
         IWF = IWC + N
         IWB = IWF + N
         IW  = IWB + N
C
         CALL MB02PD( 'Equilibrate', 'Transpose', N, N, DWORK(IU), N2,
     $                S(NP1,1), LDS, IWORK, EQUED, DWORK(IWR),
     $                DWORK(IWC), DWORK(IU+N), N2, X, LDX, RCONDU,
     $                DWORK(IWF), DWORK(IWB), IWORK(NP1), DWORK(IW),
     $                IERR )
         IF( JOBX ) THEN
C
C           Restore U(2,1) back in-situ.
C
            DO 40 I = 1, N - 1
               CALL DSWAP( N-I, DWORK(IU+N+I*(N2+1)-1), N2,
     $                     DWORK(IU+N+(I-1)*(N2+1)+1), 1 )
   40       CONTINUE
C
            IF( .NOT.LSAME( EQUED, 'N' ) ) THEN
C
C              Undo the equilibration of U(1,1) and U(2,1).
C
               ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
               COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
C
               IF( ROWEQU ) THEN
C
                  DO 60 I = 1, N
                     DWORK(IWR+I-1) = ONE / DWORK(IWR+I-1)
   60             CONTINUE
C
                  CALL MB01SD( 'Row scaling', N, N, DWORK(IU), N2,
     $                         DWORK(IWR), DWORK(IWC) )
               END IF
C
               IF( COLEQU ) THEN
C
                  DO 80 I = 1, N
                     DWORK(IWC+I-1) = ONE / DWORK(IWC+I-1)
   80             CONTINUE
C
                  CALL MB01SD( 'Column scaling', N, N, DWORK(IU), N2,
     $                         DWORK(IWR), DWORK(IWC) )
                  CALL MB01SD( 'Column scaling', N, N, DWORK(IU+N), N2,
     $                         DWORK(IWR), DWORK(IWC) )
               END IF
            END IF
C
C           Set S(2,1) to zero.
C
            CALL DLASET( 'Full', N, N, ZERO, ZERO, S(NP1,1), LDS )
         END IF
C
         PIVOTU = DWORK(IW)
C
         IF ( IERR.GT.0 ) THEN
C
C           Singular matrix. Set INFO and DWORK for error return.
C
            INFO = 5
            GO TO 160
         END IF
C
C        Make sure the solution matrix X is symmetric.
C
         DO 100 I = 1, N - 1
            CALL DAXPY( N-I, ONE, X(I,I+1), LDX, X(I+1,I), 1 )
            CALL DSCAL( N-I, HALF, X(I+1,I), 1 )
            CALL DCOPY( N-I, X(I+1,I), 1, X(I,I+1), LDX )
  100    CONTINUE
C
         IF( LSCAL ) THEN
C
C           Undo scaling for the solution matrix.
C
            IF( LSCL )
     $         CALL DLASCL( 'G', 0, 0, GNORM, QNORM, N, N, X, LDX,
     $                      IERR )
         END IF
      END IF
C
      IF ( .NOT.JOBX ) THEN
         IF ( .NOT.JOBA )
     $      WRKOPT = 0
C
C        Estimate the conditioning and compute an error bound on the
C        solution of the algebraic Riccati equation.
C
         IW = 6
         LOFACT = FACT
         IF ( NOFACT .AND. .NOT.UPDATE ) THEN
C
C           Compute Ac and its Schur factorization.
C
            IF ( DISCR ) THEN
               CALL DLASET( 'Full', N, N, ZERO, ONE, DWORK(IW), N )
               CALL DSYMM(  'Left', UPLO, N, N, ONE, G, LDG, X, LDX,
     $                      ONE, DWORK(IW), N )
               IF ( NOTRNA ) THEN
C
C                 Compute Ac = inv(I_n + G*X)*A.
C
                  CALL DLACPY( 'Full', N, N, A, LDA, T, LDT )
                  CALL DGESV( N, N, DWORK(IW), N, IWORK, T, LDT, IERR )
               ELSE
C
C                 Compute Ac = A*inv(I_n + X*G).
C
                  CALL MA02AD( 'Full', N, N, A, LDA, T, LDT )
                  CALL DGESV( N, N, DWORK(IW), N, IWORK, T, LDT, IERR )
                  DO 120 I = 2, N
                     CALL DSWAP( I-1, T(1,I), 1, T(I,1), LDT )
  120             CONTINUE
               END IF
C
            ELSE
C
               CALL DLACPY( 'Full', N, N, A, LDA, T, LDT )
               IF ( NOTRNA ) THEN
C
C                 Compute Ac = A - G*X.
C
                  CALL DSYMM( 'Left', UPLO, N, N, -ONE, G, LDG, X, LDX,
     $                        ONE, T, LDT )
               ELSE
C
C                 Compute Ac = A - X*G.
C
                  CALL DSYMM( 'Right', UPLO, N, N, -ONE, G, LDG, X, LDX,
     $                        ONE, T, LDT )
               END IF
            END IF
C
C           Compute the Schur factorization of Ac, Ac = V*T*V'.
C           Workspace:  need   5 + 5*N.
C                       prefer larger.
C
            IWR = IW
            IWI = IWR + N
            IW  = IWI + N
            LDW = LDWORK - IW + 1
C
            CALL DGEES( 'Vectors', 'Not ordered', SB02MS, N, T, LDT,
     $                  NROT, DWORK(IWR), DWORK(IWI), V, LDV, DWORK(IW),
     $                  LDW, BWORK, IERR )
C
            IF( IERR.NE.0 ) THEN
               INFO = 6
               GO TO 160
            END IF
C
            WRKOPT = MAX( WRKOPT, DWORK(IW) + DBLE( IW - 1 ) )
            LOFACT = 'F'
            IW = 6
         END IF
C
         IF ( .NOT.UPDATE ) THEN
C
C           Update G, Q, and X using the orthogonal matrix V.
C
            TRANAT = 'T'
C
C           Save the diagonal elements of G and Q.
C
            CALL DCOPY( N, G, LDG+1, DWORK(IW), 1 )
            CALL DCOPY( N, Q, LDQ+1, DWORK(IW+N), 1 )
            IW = IW + N2
C
            IF ( JOBA )
     $         CALL DLACPY( 'Full', N, N, X, LDX, S(NP1,1), LDS )
            CALL MB01RU( UPLO, TRANAT, N, N, ZERO, ONE, X, LDX, V, LDV,
     $                   X, LDX, DWORK(IW), NN, IERR )
            CALL DSCAL( N, HALF, X, LDX+1 )
            CALL MA02ED( UPLO, N, X, LDX )
            IF( .NOT.DISCR ) THEN
               CALL MA02ED( UPLO, N, G, LDG )
               CALL MA02ED( UPLO, N, Q, LDQ )
            END IF
            CALL MB01RU( UPLO, TRANAT, N, N, ZERO, ONE, G, LDG, V, LDV,
     $                   G, LDG, DWORK(IW), NN, IERR )
            CALL DSCAL( N, HALF, G, LDG+1 )
            CALL MB01RU( UPLO, TRANAT, N, N, ZERO, ONE, Q, LDQ, V, LDV,
     $                   Q, LDQ, DWORK(IW), NN, IERR )
            CALL DSCAL( N, HALF, Q, LDQ+1 )
         END IF
C
C        Estimate the conditioning and/or the error bound.
C        Workspace: 5 + MAX(1,LWS,LWE) + LWN, where
C
C           LWS = 0,       if FACT = 'F' or  LYAPUN = 'R';
C               = 5*N,     if FACT = 'N' and LYAPUN = 'O' and DICO = 'C'
C                                                         and JOB = 'C';
C               = 5*N+N*N, if FACT = 'N' and LYAPUN = 'O' and DICO = 'C'
C                                          and (JOB = 'E' or JOB = 'A');
C               = 5*N+N*N, if FACT = 'N' and LYAPUN = 'O' and
C                                                         DICO = 'D';
C           LWE = 2*N*N,                if DICO = 'C' and  JOB = 'C';
C               = 4*N*N,                if DICO = 'C' and (JOB = 'E' or
C                                                          JOB = 'A');
C               = MAX(3,2*N*N) + N*N,   if DICO = 'D' and  JOB = 'C';
C               = MAX(3,2*N*N) + 2*N*N, if DICO = 'D' and (JOB = 'E' or
C                                                          JOB = 'A');
C           LWN = 0,   if LYAPUN = 'O' or   JOB = 'C';
C               = 2*N, if LYAPUN = 'R' and DICO = 'C' and (JOB = 'E' or
C                                                          JOB = 'A');
C               = 3*N, if LYAPUN = 'R' and DICO = 'D' and (JOB = 'E' or
C                                                          JOB = 'A').
C
         LDW = LDWORK - IW + 1
         IF ( JOBA ) THEN
            JOBS = 'B'
         ELSE
            JOBS = JOB
         END IF
C
         IF ( DISCR ) THEN
            CALL SB02SD( JOBS, LOFACT, TRANA, UPLO, LYAPUN, N, A, LDA,
     $                   T, LDT, V, LDV, G, LDG, Q, LDQ, X, LDX, SEP,
     $                   RCOND, FERR, IWORK, DWORK(IW), LDW, IERR )
         ELSE
            CALL SB02QD( JOBS, LOFACT, TRANA, UPLO, LYAPUN, N, A, LDA,
     $                   T, LDT, V, LDV, G, LDG, Q, LDQ, X, LDX, SEP,
     $                   RCOND, FERR, IWORK, DWORK(IW), LDW, IERR )
         END IF
C
         WRKOPT = MAX( WRKOPT, DWORK(IW) + DBLE( IW - 1 ) )
         IF( IERR.EQ.NP1 ) THEN
            INFO = 7
         ELSE IF( IERR.GT.0 ) THEN
            INFO = 6
            GO TO 160
         END IF
C
         IF ( .NOT.UPDATE ) THEN
C
C           Restore X, G, and Q and set S(2,1) to zero, if needed.
C
            IF ( JOBA ) THEN
               CALL DLACPY( 'Full', N, N, S(NP1,1), LDS, X, LDX )
               CALL DLASET( 'Full', N, N, ZERO, ZERO, S(NP1,1), LDS )
            ELSE
               CALL MB01RU( UPLO, TRANA, N, N, ZERO, ONE, X, LDX, V,
     $                      LDV, X, LDX, DWORK(IW), NN, IERR )
               CALL DSCAL( N, HALF, X, LDX+1 )
               CALL MA02ED( UPLO, N, X, LDX )
            END IF
            IF ( LUPLO ) THEN
               LOUP = 'L'
            ELSE
               LOUP = 'U'
            END IF
C
            IW = 6
            CALL DCOPY( N, DWORK(IW), 1, G, LDG+1 )
            CALL MA02ED( LOUP, N, G, LDG )
            CALL DCOPY( N, DWORK(IW+N), 1, Q, LDQ+1 )
            CALL MA02ED( LOUP, N, Q, LDQ )
         END IF
C
      END IF
C
C     Set the optimal workspace and other details.
C
      DWORK(1) = WRKOPT
  160 CONTINUE
      IF( JBXA ) THEN
         DWORK(2) = RCONDU
         DWORK(3) = PIVOTU
         IF ( DISCR ) THEN
            DWORK(4) = RCONDA
            DWORK(5) = PIVOTA
         END IF
         IF( JOBX ) THEN
            IF ( LSCL ) THEN
               SEP = QNORM / GNORM
            ELSE
               SEP = ONE
            END IF
         END IF
      END IF
C
      RETURN
C *** Last line of SB02RD ***
      END