```--- a
+++ b/src/SB02QD.f
@@ -0,0 +1,804 @@
+      SUBROUTINE SB02QD( JOB, FACT, TRANA, UPLO, LYAPUN, N, A, LDA, T,
+     \$                   LDT, U, LDU, G, LDG, Q, LDQ, X, LDX, SEP,
+     \$                   RCOND, FERR, IWORK, DWORK, LDWORK, 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
+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
+C     PURPOSE
+C
+C     To estimate the conditioning and compute an error bound on the
+C     solution of the real continuous-time matrix algebraic Riccati
+C     equation
+C
+C         op(A)'*X + X*op(A) + Q - X*G*X = 0,                        (1)
+C
+C     where op(A) = A or A' (A**T) and Q, G are symmetric (Q = Q**T,
+C     G = G**T). The matrices A, Q and G are N-by-N and the solution X
+C     is N-by-N.
+C
+C     ARGUMENTS
+C
+C     Mode Parameters
+C
+C     JOB     CHARACTER*1
+C             Specifies the computation to be performed, as follows:
+C             = 'C':  Compute the reciprocal condition number only;
+C             = 'E':  Compute the error bound only;
+C             = 'B':  Compute both the reciprocal condition number and
+C                     the error bound.
+C
+C     FACT    CHARACTER*1
+C             Specifies whether or not the real Schur factorization of
+C             the matrix Ac = A - G*X (if TRANA = 'N') or Ac = A - X*G
+C             (if TRANA = 'T' or 'C') is supplied on entry, as follows:
+C             = 'F':  On entry, T and U (if LYAPUN = 'O') contain the
+C                     factors from the real Schur factorization of the
+C                     matrix Ac;
+C             = 'N':  The Schur factorization of Ac will be computed
+C                     and the factors will be stored in T and U (if
+C                     LYAPUN = 'O').
+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 part of the symmetric matrices Q and G is
+C             to be used, as follows:
+C             = 'U':  Upper triangular part;
+C             = 'L':  Lower triangular part.
+C
+C     LYAPUN  CHARACTER*1
+C             Specifies whether or not the original Lyapunov equations
+C             should be solved in the iterative estimation process,
+C             as follows:
+C             = 'O':  Solve the original Lyapunov equations, updating
+C                     the right-hand sides and solutions with the
+C                     matrix U, e.g., RHS <-- U'*RHS*U;
+C             = 'R':  Solve reduced Lyapunov equations only, without
+C                     updating the right-hand sides and solutions.
+C
+C     Input/Output Parameters
+C
+C     N       (input) INTEGER
+C             The order of the matrices A, X, Q, and G.  N >= 0.
+C
+C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
+C             If FACT = 'N' or LYAPUN = 'O', the leading N-by-N part of
+C             this array must contain the matrix A.
+C             If FACT = 'F' and LYAPUN = 'R', A is not referenced.
+C
+C     LDA     INTEGER
+C             The leading dimension of the array A.
+C             LDA >= max(1,N), if FACT = 'N' or  LYAPUN = 'O';
+C             LDA >= 1,        if FACT = 'F' and LYAPUN = 'R'.
+C
+C     T       (input or output) DOUBLE PRECISION array, dimension
+C             (LDT,N)
+C             If FACT = 'F', then T is an input argument and on entry,
+C             the leading N-by-N upper Hessenberg part of this array
+C             must contain the upper quasi-triangular matrix T in Schur
+C             canonical form from a Schur factorization of Ac (see
+C             argument FACT).
+C             If FACT = 'N', then T is an output argument and on exit,
+C             if INFO = 0 or INFO = N+1, the leading N-by-N upper
+C             Hessenberg part of this array contains the upper quasi-
+C             triangular matrix T in Schur canonical form from a Schur
+C             factorization of Ac (see argument FACT).
+C
+C     LDT     INTEGER
+C             The leading dimension of the array T.  LDT >= max(1,N).
+C
+C     U       (input or output) DOUBLE PRECISION array, dimension
+C             (LDU,N)
+C             If LYAPUN = 'O' and FACT = 'F', then U is an input
+C             argument and on entry, the leading N-by-N part of this
+C             array must contain the orthogonal matrix U from a real
+C             Schur factorization of Ac (see argument FACT).
+C             If LYAPUN = 'O' and FACT = 'N', then U is an output
+C             argument and on exit, if INFO = 0 or INFO = N+1, it
+C             contains the orthogonal N-by-N matrix from a real Schur
+C             factorization of Ac (see argument FACT).
+C             If LYAPUN = 'R', the array U is not referenced.
+C
+C     LDU     INTEGER
+C             The leading dimension of the array U.
+C             LDU >= 1,        if LYAPUN = 'R';
+C             LDU >= MAX(1,N), if LYAPUN = 'O'.
+C
+C     G       (input) DOUBLE PRECISION array, dimension (LDG,N)
+C             If UPLO = 'U', the leading N-by-N upper triangular part of
+C             this array must contain the upper triangular part of the
+C             matrix G.
+C             If UPLO = 'L', the leading N-by-N lower triangular part of
+C             this array must contain the lower triangular part of the
+C             matrix G.                     _
+C             Matrix G should correspond to G in the "reduced" Riccati
+C             equation (with matrix T, instead of A), if LYAPUN = 'R'.
+C             See METHOD.
+C
+C     LDG     INTEGER
+C             The leading dimension of the array G.  LDG >= max(1,N).
+C
+C     Q       (input) DOUBLE PRECISION array, dimension (LDQ,N)
+C             If UPLO = 'U', the leading N-by-N upper triangular part of
+C             this array must contain the upper triangular part of the
+C             matrix Q.
+C             If UPLO = 'L', the leading N-by-N lower triangular part of
+C             this array must contain the lower triangular part of the
+C             matrix Q.                     _
+C             Matrix Q should correspond to Q in the "reduced" Riccati
+C             equation (with matrix T, instead of A), if LYAPUN = 'R'.
+C             See METHOD.
+C
+C     LDQ     INTEGER
+C             The leading dimension of the array Q.  LDQ >= max(1,N).
+C
+C     X       (input) DOUBLE PRECISION array, dimension (LDX,N)
+C             The leading N-by-N part of this array must contain the
+C             symmetric solution matrix of the original Riccati
+C             equation (with matrix A), if LYAPUN = 'O', or of the
+C             "reduced" Riccati equation (with matrix T), if
+C             LYAPUN = 'R'. See METHOD.
+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 = 'B', the estimated quantity
+C             sep(op(Ac),-op(Ac)').
+C             If N = 0, or X = 0, or JOB = 'E', SEP is not referenced.
+C
+C     RCOND   (output) DOUBLE PRECISION
+C             If JOB = 'C' or JOB = 'B', an estimate of the reciprocal
+C             condition number of the continuous-time Riccati equation.
+C             If N = 0 or X = 0, RCOND is set to 1 or 0, respectively.
+C             If JOB = 'E', RCOND is not referenced.
+C
+C     FERR    (output) DOUBLE PRECISION
+C             If JOB = 'E' or JOB = 'B', an estimated forward error
+C             bound for the solution X. If XTRUE is the true solution,
+C             FERR bounds the magnitude of the largest entry in
+C             (X - XTRUE) divided by the magnitude of the largest entry
+C             in X.
+C             If N = 0 or X = 0, FERR is set to 0.
+C             If JOB = 'C', FERR is not referenced.
+C
+C     Workspace
+C
+C     IWORK   INTEGER array, dimension (N*N)
+C
+C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
+C             On exit, if INFO = 0 or INFO = N+1, DWORK(1) returns the
+C             optimal value of LDWORK.
+C
+C     LDWORK  INTEGER
+C             The dimension of the array DWORK.
+C             Let LWA = N*N, if LYAPUN = 'O' and JOB = 'E' or 'B';
+C                 LWA = 0,   otherwise.
+C             If FACT = 'N', then
+C                LDWORK  = MAX(1, 5*N, 2*N*N),        if JOB = 'C';
+C                LDWORK  = MAX(1, LWA + 5*N, 4*N*N ), if JOB = 'E', 'B'.
+C             If FACT = 'F', then
+C                LDWORK  = MAX(1, 2*N*N),  if JOB = 'C';
+C                LDWORK  = MAX(1, 4*N*N ), if JOB = 'E' or 'B'.
+C             For good performance, LDWORK must generally be larger.
+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             > 0:  if INFO = i, i <= N, the QR algorithm failed to
+C                   complete the reduction of the matrix Ac to Schur
+C                   canonical form (see LAPACK Library routine DGEES);
+C                   on exit, the matrix T(i+1:N,i+1:N) contains the
+C                   partially converged Schur form, and DWORK(i+1:N) and
+C                   DWORK(N+i+1:2*N) contain the real and imaginary
+C                   parts, respectively, of the converged eigenvalues;
+C                   this error is unlikely to appear;
+C             = N+1:  if the matrices T and -T' have common or very
+C                   close eigenvalues; perturbed values were used to
+C                   solve Lyapunov equations, but the matrix T, if given
+C                   (for FACT = 'F'), is unchanged.
+C
+C     METHOD
+C
+C     The condition number of the 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     and Ac = A - G*X (if TRANA = 'N') or Ac = A - X*G (if TRANA = 'T'
+C     or 'C'). Note that the Riccati equation (1) is equivalent to
+C                _   _         _   _ _ _
+C         op(T)'*X + X*op(T) + Q + X*G*X = 0,                        (2)
+C           _           _               _
+C     where X = U'*X*U, Q = U'*Q*U, and G = U'*G*U, with U the
+C     orthogonal matrix reducing Ac to a real Schur form, T = U'*Ac*U.
+C
+C     The routine estimates the quantities
+C
+C     sep(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 [2].
+C
+C     REFERENCES
+C
+C     [1] 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     [2] Higham, N.J.
+C         Perturbation theory and backward error for AX-XB=C.
+C         BIT, vol. 33, pp. 124-136, 1993.
+C
+C     [3] 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.
+C     The accuracy of the estimates obtained depends on the solution
+C     accuracy and on the properties of the 1-norm estimator.
+C
+C
+C     The option LYAPUN = 'R' may occasionally produce slightly worse
+C     or better estimates, and it is much faster than the option 'O'.
+C     When SEP is computed and it is zero, the routine returns
+C     immediately, with RCOND and FERR (if requested) set to 0 and 1,
+C     respectively. In this case, the equation is singular.
+C
+C     CONTRIBUTOR
+C
+C     P.Hr. Petkov, Technical University of Sofia, December 1998.
+C     V. Sima, Katholieke Univ. Leuven, Belgium, February 1999.
+C
+C     REVISIONS
+C
+C     V. Sima, Research Institute for Informatics, Bucharest, Oct. 2004.
+C
+C     KEYWORDS
+C
+C     Conditioning, error estimates, orthogonal transformation,
+C     real Schur form, Riccati equation.
+C
+C     ******************************************************************
+C
+C     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE, TWO, FOUR, HALF
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
+     \$                     FOUR = 4.0D+0, HALF = 0.5D+0 )
+C     ..
+C     .. Scalar Arguments ..
+      CHARACTER          FACT, JOB, LYAPUN, TRANA, UPLO
+      INTEGER            INFO, LDA, LDG, LDQ, LDT, LDU, LDWORK, LDX, N
+      DOUBLE PRECISION   FERR, RCOND, SEP
+C     ..
+C     .. Array Arguments ..
+      INTEGER            IWORK( * )
+      DOUBLE PRECISION   A( LDA, * ), DWORK( * ),  G( LDG, * ),
+     \$                   Q( LDQ, * ), T( LDT, * ), U( LDU, * ),
+     \$                   X( LDX, * )
+C     ..
+C     .. Local Scalars ..
+      LOGICAL            JOBB, JOBC, JOBE, LOWER, NEEDAC, NOFACT,
+     \$                   NOTRNA, UPDATE
+      CHARACTER          LOUP, SJOB, TRANAT
+      INTEGER            I, IABS, INFO2, IRES, ITMP, IXBS, J, JJ, JX,
+     \$                   KASE, LDW, LWA, NN, SDIM, WRKOPT
+      DOUBLE PRECISION   ANORM, BIGNUM, DENOM, EPS, EPSN, EST, GNORM,
+     \$                   PINORM, QNORM, SCALE, SIG, TEMP, THNORM, TMAX,
+     \$                   XANORM, XNORM
+C     ..
+C     .. Local Arrays ..
+      LOGICAL            BWORK( 1 )
+C     ..
+C     .. External Functions ..
+      LOGICAL            LSAME, SELECT
+      DOUBLE PRECISION   DLAMCH, DLANGE, DLANHS, DLANSY
+      EXTERNAL           DLAMCH, DLANGE, DLANHS, DLANSY, LSAME, SELECT
+C     ..
+C     .. External Subroutines ..
+      EXTERNAL           DAXPY, DCOPY, DGEES, DLACON, DLACPY, DSCAL,
+     \$                   DSYMM, DSYR2K, MA02ED, MB01RU, MB01UD, SB03MY,
+     \$                   SB03QX, SB03QY, XERBLA
+C     ..
+C     .. Intrinsic Functions ..
+      INTRINSIC          ABS, DBLE, INT, MAX, MIN
+C     ..
+C     .. Executable Statements ..
+C
+C     Decode and Test input parameters.
+C
+      JOBC   = LSAME( JOB,    'C' )
+      JOBE   = LSAME( JOB,    'E' )
+      JOBB   = LSAME( JOB,    'B' )
+      NOFACT = LSAME( FACT,   'N' )
+      NOTRNA = LSAME( TRANA,  'N' )
+      LOWER  = LSAME( UPLO,   'L' )
+      UPDATE = LSAME( LYAPUN, 'O' )
+C
+      NEEDAC = UPDATE .AND. .NOT.JOBC
+C
+      NN = N*N
+      IF( NEEDAC ) THEN
+         LWA = NN
+      ELSE
+         LWA = 0
+      END IF
+C
+      IF( NOFACT ) THEN
+         IF( JOBC ) THEN
+            LDW = MAX( 5*N, 2*NN )
+         ELSE
+            LDW = MAX( LWA + 5*N, 4*NN )
+         END IF
+      ELSE
+         IF( JOBC ) THEN
+            LDW = 2*NN
+         ELSE
+            LDW = 4*NN
+         END IF
+      END IF
+C
+      INFO = 0
+      IF( .NOT.( JOBB .OR. JOBC .OR. JOBE ) ) THEN
+         INFO = -1
+      ELSE IF( .NOT.( NOFACT .OR. LSAME( FACT,   'F' ) ) ) THEN
+         INFO = -2
+      ELSE IF( .NOT.( NOTRNA .OR. LSAME( TRANA,  'T' ) .OR.
+     \$                            LSAME( TRANA,  'C' ) ) ) THEN
+         INFO = -3
+      ELSE IF( .NOT.( LOWER  .OR. LSAME( UPLO,   'U' ) ) ) THEN
+         INFO = -4
+      ELSE IF( .NOT.( UPDATE .OR. LSAME( LYAPUN, 'R' ) ) ) THEN
+         INFO = -5
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -6
+      ELSE IF( LDA.LT.1 .OR.
+     \$       ( LDA.LT.N .AND. ( UPDATE .OR. NOFACT ) ) ) THEN
+         INFO = -8
+      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+         INFO = -10
+      ELSE IF( LDU.LT.1 .OR. ( LDU.LT.N .AND. UPDATE ) ) THEN
+         INFO = -12
+      ELSE IF( LDG.LT.MAX( 1, N ) ) THEN
+         INFO = -14
+      ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
+         INFO = -16
+      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
+         INFO = -18
+      ELSE IF( LDWORK.LT.MAX( 1, LDW ) ) THEN
+         INFO = -24
+      END IF
+C
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SB02QD', -INFO )
+         RETURN
+      END IF
+C
+C     Quick return if possible.
+C
+      IF( N.EQ.0 ) THEN
+         IF( .NOT.JOBE )
+     \$      RCOND = ONE
+         IF( .NOT.JOBC )
+     \$      FERR  = ZERO
+         DWORK( 1 ) = ONE
+         RETURN
+      END IF
+C
+C     Compute the 1-norm of the matrix X.
+C
+      XNORM = DLANSY( '1-norm', UPLO, N, X, LDX, DWORK )
+      IF( XNORM.EQ.ZERO ) THEN
+C
+C        The solution is zero.
+C
+         IF( .NOT.JOBE )
+     \$      RCOND = ZERO
+         IF( .NOT.JOBC )
+     \$      FERR  = ZERO
+         DWORK( 1 ) = DBLE( N )
+         RETURN
+      END IF
+C
+C     Workspace usage.
+C
+      IXBS = 0
+      ITMP = IXBS + NN
+      IABS = ITMP + NN
+      IRES = IABS + NN
+C
+C     Workspace:  LWR, where
+C                 LWR = N*N, if LYAPUN = 'O' and JOB = 'E' or 'B', or
+C                               FACT = 'N',
+C                 LWR = 0,   otherwise.
+C
+      IF( NEEDAC .OR. NOFACT ) THEN
+C
+         CALL DLACPY( 'Full', N, N, A, LDA, DWORK, N )
+         IF( NOTRNA ) THEN
+C
+C           Compute Ac = A - G*X.
+C
+            CALL DSYMM( 'Left', UPLO, N, N, -ONE, G, LDG, X, LDX, ONE,
+     \$                  DWORK, N )
+         ELSE
+C
+C           Compute Ac = A - X*G.
+C
+            CALL DSYMM( 'Right', UPLO, N, N, -ONE, G, LDG, X, LDX, ONE,
+     \$                  DWORK, N )
+         END IF
+C
+         WRKOPT = DBLE( NN )
+         IF( NOFACT )
+     \$      CALL DLACPY( 'Full', N, N, DWORK, N, T, LDT )
+      ELSE
+         WRKOPT = DBLE( N )
+      END IF
+C
+      IF( NOFACT ) THEN
+C
+C        Compute the Schur factorization of Ac, Ac = U*T*U'.
+C        Workspace:  need   LWA + 5*N;
+C                    prefer larger;
+C                    LWA = N*N, if LYAPUN = 'O' and JOB = 'E' or 'B';
+C                    LWA = 0,   otherwise.
+C        (Note: Comments in the code beginning "Workspace:" describe the
+C        minimal amount of real workspace needed at that point in the
+C        code, as well as the preferred amount for good performance.)
+C
+         IF( UPDATE ) THEN
+            SJOB = 'V'
+         ELSE
+            SJOB = 'N'
+         END IF
+         CALL DGEES( SJOB, 'Not ordered', SELECT, N, T, LDT, SDIM,
+     \$               DWORK( LWA+1 ), DWORK( LWA+N+1 ), U, LDU,
+     \$               DWORK( LWA+2*N+1 ), LDWORK-LWA-2*N, BWORK, INFO )
+         IF( INFO.GT.0 ) THEN
+            IF( LWA.GT.0 )
+     \$         CALL DCOPY( 2*N, DWORK( LWA+1 ), 1, DWORK, 1 )
+            RETURN
+         END IF
+C
+         WRKOPT = MAX( WRKOPT, INT( DWORK( LWA+2*N+1 ) ) + LWA + 2*N )
+      END IF
+      IF( NEEDAC )
+     \$   CALL DLACPY( 'Full', N, N, DWORK, N, DWORK( IABS+1 ), N )
+C
+      IF( NOTRNA ) THEN
+         TRANAT = 'T'
+      ELSE
+         TRANAT = 'N'
+      END IF
+C
+      IF( .NOT.JOBE ) THEN
+C
+C        Estimate sep(op(Ac),-op(Ac)') = sep(op(T),-op(T)') and
+C        norm(Theta).
+C        Workspace LWA + 2*N*N.
+C
+         CALL SB03QY( 'Both', TRANA, LYAPUN, N, T, LDT, U, LDU, X, LDX,
+     \$                SEP, THNORM, IWORK, DWORK, LDWORK, INFO )
+C
+         WRKOPT = MAX( WRKOPT, LWA + 2*NN )
+C
+C        Return if the equation is singular.
+C
+         IF( SEP.EQ.ZERO ) THEN
+            RCOND = ZERO
+            IF( JOBB )
+     \$         FERR = ONE
+            DWORK( 1 ) = DBLE( WRKOPT )
+            RETURN
+         END IF
+C
+C        Estimate norm(Pi).
+C        Workspace LWA + 2*N*N.
+C
+         KASE = 0
+C
+C        REPEAT
+   10    CONTINUE
+         CALL DLACON( NN, DWORK( ITMP+1 ), DWORK, IWORK, EST, KASE )
+         IF( KASE.NE.0 ) THEN
+C
+C           Select the triangular part of symmetric matrix to be used.
+C
+            IF( DLANSY( '1-norm', 'Upper', N, DWORK, N, DWORK( ITMP+1 ))
+     \$          .GE.
+     \$          DLANSY( '1-norm', 'Lower', N, DWORK, N, DWORK( ITMP+1 ))
+     \$        ) THEN
+               LOUP = 'U'
+            ELSE
+               LOUP = 'L'
+            END IF
+C
+C           Compute RHS = X*W*X.
+C
+            CALL MB01RU( LOUP, 'No Transpose', N, N, ZERO, ONE, DWORK,
+     \$                   N, X, LDX, DWORK, N, DWORK( ITMP+1 ), NN,
+     \$                   INFO2 )
+            CALL DSCAL( N, HALF, DWORK, N+1 )
+C
+            IF( UPDATE ) THEN
+C
+C              Transform the right-hand side: RHS := U'*RHS*U.
+C
+               CALL MB01RU( LOUP, 'Transpose', N, N, ZERO, ONE, DWORK,
+     \$                      N, U, LDU, DWORK, N, DWORK( ITMP+1 ), NN,
+     \$                      INFO2 )
+               CALL DSCAL( N, HALF, DWORK, N+1 )
+            END IF
+C
+C           Fill in the remaining triangle of the symmetric matrix.
+C
+            CALL MA02ED( LOUP, N, DWORK, N )
+C
+            IF( KASE.EQ.1 ) THEN
+C
+C              Solve op(T)'*Y + Y*op(T) = scale*RHS.
+C
+               CALL SB03MY( TRANA, N, T, LDT, DWORK, N, SCALE, INFO2 )
+            ELSE
+C
+C              Solve op(T)*W + W*op(T)' = scale*RHS.
+C
+               CALL SB03MY( TRANAT, N, T, LDT, DWORK, N, SCALE, INFO2 )
+            END IF
+C
+            IF( UPDATE ) THEN
+C
+C              Transform back to obtain the solution: Z := U*Z*U', with
+C              Z = Y or Z = W.
+C
+               CALL MB01RU( LOUP, 'No transpose', N, N, ZERO, ONE,
+     \$                      DWORK, N, U, LDU, DWORK, N, DWORK( ITMP+1 ),
+     \$                      NN, INFO2 )
+               CALL DSCAL( N, HALF, DWORK, N+1 )
+C
+C              Fill in the remaining triangle of the symmetric matrix.
+C
+               CALL MA02ED( LOUP, N, DWORK, N )
+            END IF
+            GO TO 10
+         END IF
+C        UNTIL KASE = 0
+C
+         IF( EST.LT.SCALE ) THEN
+            PINORM = EST / SCALE
+         ELSE
+            BIGNUM = ONE / DLAMCH( 'Safe minimum' )
+            IF( EST.LT.SCALE*BIGNUM ) THEN
+               PINORM = EST / SCALE
+            ELSE
+               PINORM = BIGNUM
+            END IF
+         END IF
+C
+C        Compute the 1-norm of A or T.
+C
+         IF( UPDATE ) THEN
+            ANORM = DLANGE( '1-norm', N, N, A, LDA, DWORK )
+         ELSE
+            ANORM = DLANHS( '1-norm', N, T, LDT, DWORK )
+         END IF
+C
+C        Compute the 1-norms of the matrices Q and G.
+C
+         QNORM = DLANSY( '1-norm', UPLO, N, Q, LDQ, DWORK )
+         GNORM = DLANSY( '1-norm', UPLO, N, G, LDG, DWORK )
+C
+C        Estimate the reciprocal condition number.
+C
+         TMAX = MAX( SEP, XNORM, ANORM, GNORM )
+         IF( TMAX.LE.ONE ) THEN
+            TEMP  = SEP*XNORM
+            DENOM = QNORM + ( SEP*ANORM )*THNORM +
+     \$                      ( SEP*GNORM )*PINORM
+         ELSE
+            TEMP  =   ( SEP / TMAX )*( XNORM / TMAX )
+            DENOM = ( ( ONE / TMAX )*( QNORM / TMAX ) ) +
+     \$              ( ( SEP / TMAX )*( ANORM / TMAX ) )*THNORM +
+     \$              ( ( SEP / TMAX )*( GNORM / TMAX ) )*PINORM
+         END IF
+         IF( TEMP.GE.DENOM ) THEN
+            RCOND = ONE
+         ELSE
+            RCOND = TEMP / DENOM
+         END IF
+      END IF
+C
+      IF( .NOT.JOBC ) THEN
+C
+C        Form a triangle of the residual matrix
+C          R = op(A)'*X + X*op(A) + Q - X*G*X,
+C        or           _   _         _   _ _ _
+C          R = op(T)'*X + X*op(T) + Q + X*G*X,
+C        exploiting the symmetry.
+C        Workspace 4*N*N.
+C
+         IF( UPDATE ) THEN
+            CALL DLACPY( UPLO, N, N, Q, LDQ, DWORK( IRES+1 ), N )
+            CALL DSYR2K( UPLO, TRANAT, N, N, ONE, A, LDA, X, LDX, ONE,
+     \$                   DWORK( IRES+1 ), N )
+            SIG = -ONE
+         ELSE
+            CALL MB01UD( 'Right', TRANA, N, N, ONE, T, LDT, X, LDX,
+     \$                   DWORK( IRES+1 ), N, INFO2 )
+            JJ = IRES + 1
+            IF( LOWER ) THEN
+               DO 20 J = 1, N
+                  CALL DAXPY( N-J+1, ONE, DWORK( JJ ), N, DWORK( JJ ),
+     \$                        1 )
+                  CALL DAXPY( N-J+1, ONE, Q( J, J ), 1, DWORK( JJ ), 1 )
+                  JJ = JJ + N + 1
+   20          CONTINUE
+            ELSE
+               DO 30 J = 1, N
+                  CALL DAXPY( J, ONE, DWORK( IRES+J ), N, DWORK( JJ ),
+     \$                        1 )
+                  CALL DAXPY( J, ONE, Q( 1, J ), 1, DWORK( JJ ), 1 )
+                  JJ = JJ + N
+   30          CONTINUE
+            END IF
+            SIG = ONE
+         END IF
+         CALL MB01RU( UPLO, TRANAT, N, N, ONE, SIG, DWORK( IRES+1 ),
+     \$                N, X, LDX, G, LDG, DWORK( ITMP+1 ), NN, INFO2 )
+C
+C        Get the machine precision.
+C
+         EPS  = DLAMCH( 'Epsilon' )
+         EPSN = EPS*DBLE( N + 4 )
+         TEMP = EPS*FOUR
+C
+C        Add to abs(R) a term that takes account of rounding errors in
+C        forming R:
+C         abs(R) := abs(R) + EPS*(4*abs(Q) + (n+4)*(abs(op(Ac))'*abs(X)
+C                 + abs(X)*abs(op(Ac))) + 2*(n+1)*abs(X)*abs(G)*abs(X)),
+C        or                             _                           _
+C         abs(R) := abs(R) + EPS*(4*abs(Q) + (n+4)*(abs(op(T))'*abs(X)
+C                       _                            _      _      _
+C                 + abs(X)*abs(op(T))) + 2*(n+1)*abs(X)*abs(G)*abs(X)),
+C        where EPS is the machine precision.
+C
+         DO 50 J = 1, N
+            DO 40 I = 1, N
+               DWORK( IXBS+(J-1)*N+I ) = ABS( X( I, J ) )
+   40       CONTINUE
+   50    CONTINUE
+C
+         IF( LOWER ) THEN
+            DO 70 J = 1, N
+               DO 60 I = J, N
+                  DWORK( IRES+(J-1)*N+I ) = TEMP*ABS( Q( I, J ) ) +
+     \$                   ABS( DWORK( IRES+(J-1)*N+I ) )
+   60          CONTINUE
+   70       CONTINUE
+         ELSE
+            DO 90 J = 1, N
+               DO 80 I = 1, J
+                  DWORK( IRES+(J-1)*N+I ) = TEMP*ABS( Q( I, J ) ) +
+     \$                   ABS( DWORK( IRES+(J-1)*N+I ) )
+   80          CONTINUE
+   90       CONTINUE
+         END IF
+C
+         IF( UPDATE ) THEN
+C
+            DO 110 J = 1, N
+               DO 100 I = 1, N
+                  DWORK( IABS+(J-1)*N+I ) =
+     \$               ABS( DWORK( IABS+(J-1)*N+I ) )
+  100          CONTINUE
+  110       CONTINUE
+C
+            CALL DSYR2K( UPLO, TRANAT, N, N, EPSN, DWORK( IABS+1 ), N,
+     \$                   DWORK( IXBS+1 ), N, ONE,  DWORK( IRES+1 ), N )
+         ELSE
+C
+            DO 130 J = 1, N
+               DO 120 I = 1, MIN( J+1, N )
+                  DWORK( IABS+(J-1)*N+I ) = ABS( T( I, J ) )
+  120          CONTINUE
+  130       CONTINUE
+C
+            CALL MB01UD( 'Left', TRANAT, N, N, EPSN, DWORK( IABS+1 ), N,
+     \$                   DWORK( IXBS+1), N, DWORK( ITMP+1 ), N, INFO2 )
+            JJ = IRES + 1
+            JX = ITMP + 1
+            IF( LOWER ) THEN
+               DO 140 J = 1, N
+                  CALL DAXPY( N-J+1, ONE, DWORK( JX ), N, DWORK( JX ),
+     \$                        1 )
+                  CALL DAXPY( N-J+1, ONE, DWORK( JX ), 1, DWORK( JJ ),
+     \$                        1 )
+                  JJ = JJ + N + 1
+                  JX = JX + N + 1
+  140          CONTINUE
+            ELSE
+               DO 150 J = 1, N
+                  CALL DAXPY( J, ONE, DWORK( ITMP+J ), N, DWORK( JX ),
+     \$                        1 )
+                  CALL DAXPY( J, ONE, DWORK( JX ), 1, DWORK( JJ ), 1 )
+                  JJ = JJ + N
+                  JX = JX + N
+  150          CONTINUE
+            END IF
+         END IF
+C
+         IF( LOWER ) THEN
+            DO 170 J = 1, N
+               DO 160 I = J, N
+                  DWORK( IABS+(J-1)*N+I ) = ABS( G( I, J ) )
+  160          CONTINUE
+  170       CONTINUE
+         ELSE
+            DO 190 J = 1, N
+               DO 180 I = 1, J
+                  DWORK( IABS+(J-1)*N+I ) = ABS( G( I, J ) )
+  180          CONTINUE
+  190       CONTINUE
+         END IF
+C
+         CALL MB01RU( UPLO, TRANA, N, N, ONE, EPS*DBLE( 2*( N + 1 ) ),
+     \$                DWORK( IRES+1 ), N, DWORK( IXBS+1), N,
+     \$                DWORK( IABS+1 ), N, DWORK( ITMP+1 ), NN, INFO2 )
+C
+         WRKOPT = MAX( WRKOPT, 4*NN )
+C
+C        Compute forward error bound, using matrix norm estimator.
+C        Workspace 4*N*N.
+C
+         XANORM = DLANSY( 'Max', UPLO, N, X, LDX, DWORK )
+C
+         CALL SB03QX( TRANA, UPLO, LYAPUN, N, XANORM, T, LDT, U, LDU,
+     \$                DWORK( IRES+1 ), N, FERR, IWORK, DWORK, IRES,
+     \$                INFO )
+      END IF
+C
+      DWORK( 1 ) = DBLE( WRKOPT )
+      RETURN
+C
+C *** Last line of SB02QD ***
+      END
```