# SLICOT Library Release Notes

## Version v5.9.1

Version 5.9.1 is the first version maintained by M. Köhler and J. Saak. It does
not add any new routines to SLICOT.

### Improvements

* The build system now uses CMAKE. The old makefile based system works in parallel
but will be removed in the future with the next minor release.
* CI Runs on Github and Gitlab - See Github #2
* RPM packaging, see INSTALL.md for details

### Bugfixes

* Fix: implicit variables are now declared
* Fix: stack overflow in the test suite on Windows
* Fix: Missing LSAME declaration in TMB03LD.f - See Github PR #14
* Fix: Invalid work space computation in IB01AD - See Github PR #18
* Fix: Include deprecated LAPACK routines - See Github #24 and #10


## Version v5.9

Version v5.9 of the **SLICOT** Library includes 11 new routines and 9 updated routines, in comparison to Version v5.8 of the library.

**New Routines**

**`AB13HD:`**  Computes the L_infinity-norm of a standard or descriptor system using structure-preserving eigenvalue computations.

**`MA01DD:`**  Computes an approximate symmetric chordal metric,

               d = min{ |a - b|, |1/a - 1/b| }

for two complex numbers `a` and `b`, given by real and imaginary parts.

**`MA01DZ:`**  Computes an approximate symmetric chordal metric `d` for two, possibly infinite, complex numbers `a` and `b`, given by fractions with non-negative denominators.

**`MA02RD:`**  Sorts the elements of a vector in increasing or decreasing order, and rearranges the elements of another vector using the same permutations.

**`MA02SD:`**  Computes the smallest nonzero absolute value of the elements of a real matrix.

**`MB04RD:`**  Reduces a real matrix pair `(A,B)` in generalized real Schur form to a block-diagonal form using well-conditioned non-orthogonal equivalence transformations. The condition numbers of the transformations used for reduction are roughly bounded by `PMAX`, where `PMAX >= 1` is a given real value. The transformations are optionally postmultiplied in two given matrices `X` and `Y`. The generalized Schur form is optionally ordered, so that clustered eigenvalues are grouped in the same block.

**`MB04RS:`**  Attempts to solve the generalized real Sylvester equation

               A * R - L * B = scale * C,                            (1)
               D * R - L * E = scale * F,

using Level 1 and Level 2 BLAS, but aborts the calculations when the absolute value of an element of `R` or `L` is greater than a given real number `PMAX >= 1`.

**`MB04RT:`**  Attempts to solve the generalized real Sylvester equation (1) using Level 3 BLAS, but aborts the calculations when the absolute value of an element of `R` or `L` is greater than a given real number `PMAX >= 1`.

**`MB04RV:`**  Attempts to solve the generalized complex Sylvester equation (1) using Level 1 and Level 2 BLAS, but aborts the calculations when the absolute value of an element of `R` or `L` is greater than a given real number `PMAX >= 1`.

**`MB04RW:`**  Attempts to solve the generalized complex Sylvester equation (1) using Level 3 BLAS, but aborts the calculations when the absolute value of an element of `R` or `L` is greater than a given real number `PMAX >= 1`.

**`MB04RZ:`**  Reduces a complex matrix pair `(A,B)` in generalized complex Schur form to a block-diagonal form using well-conditioned non-unitary equivalence transformations. The condition numbers of the transformations used for reduction are roughly bounded by `PMAX`, where `PMAX >= 1` is a given real value. The transformations are optionally postmultiplied in two given matrices `X` and `Y`. The generalized Schur form is optionally ordered, so that clustered eigenvalues are grouped in the same block.

**Updated Routines**

**`AB13DX:`**  Replaced `CWORK` and `LCWORK` by `ZWORK` and `LZWORK`, respectively, to agree to the **SLICOT** standards. Made several cosmetical changes.

**`AB13ID:`**  Placed the lines defining the external subroutines after the lines defining the external functions, to agree to the **SLICOT** standard. Replaced `TOLDEF` by `ZERO` in the `MB03OD` call with `LDWORK = -1`. Changed the tests for checking if the estimated rank might be incorrect. The test for setting the value 1 to `IWARN` has been modified for better significance and `IWARN` parameter description has been updated. The submatrix `R22` has also been set to zero. Used a safeguard against the huge workspace size returned by some implementations of the LAPACK routine `ZGESVD`. Made several cosmetical changes.

**`MA02JD:`**  Placed the lines defining the external subroutines after the lines defining the external functions, to agree to the **SLICOT** standard.

**`MA02JZ:`**  Placed the lines defining the external subroutines after the lines defining the external functions, to agree to the **SLICOT** standard.

**`MB02SZ:`**  Replaced the references to the `DCABS1` routine, which is not available in some `BLAS` implementations, to equivalent references to a locally defined statement functions.

**`MB4DLZ:`**  Replaced the references to the `DCABS1` routine, which is not available in some `BLAS` implementations, to equivalent references to a locally defined statement functions.

**`SB03OD:`**  Inserted a call to `DGERQF` for computing the optimal workspace when `TRANS = 'T'`.  Without this call, an error could appear when `M > N`.

**`SB03OZ:`**  Inserted a call to `ZGERQF` for computing the optimal workspace when `TRANS = 'C'`.  Without this call, an error could appear when `M > N`.

**`TG01JD:`**  Replaced the argument `ZERO` of the `TG01AD` call in the line 385 by a value computed based on the norms of the system matrices.

**Updated makefile**

The files **`makefile`** and **`makefile_Unix`** in the **`src`** sub-directory have been updated to include the new routines.

**New and Updated Documentation Files**

Documentation files have been included for the new routines, and the existing documentation files have been updated when needed.

## Version v5.8, Update 1

**Updated Routines**

**`MB01PD:`** Removed `SAVE` statement and variable `FIRST` to get a thread safe version.

**`MB03VY:`**  Set `A( ILO, ILO, J )` to 1 if `IHI = ILO` in the loop labelled 20.

**`MB04BD`**:  Changed some part of finding the number of infinite eigenvalues (replaced `DE(1,2)` by `DE(1,3)` in the line 495, and replaced `NINF = MAX( I, J )` by `NINF = MAX( I, J )/2` in the line 513).

**`MB04BP`**:  Changed some part of finding the number of infinite eigenvalues (replaced `DE(1,2)` by `DE(1,3)` in the line 509, and replaced `NINF = MAX( I, J )` by `NINF = MAX( I, J )/2` in the line 527).  Called **`MB04BD`** internally, if on entry `INFO = 0` or if `INFO < 0` and `M` is sufficiently small (`M <= NX`, with parameter `NX` set to 250).  Updated the documentation file, **`MB04BP.html`**.

**Updated makefile**

The files **`makefile`** in the **`examples`** and **`src`** sub-directories have been updated, by alphabetically ordering the names, and adding the missing ones.

**New files**

Added **`make.inc`** and **`makefile`** in the **`slicot`** directory for Windows platform with Intel Fortran compilers.

Added **`make_Unix.inc`** in **`slicot`** directory and **`makefile_Unix`** in **`slicot`** directory and subdirectories **`examples`** and **`src`**, for Unix-like platfoms with gfortran compiler.

Added the file **`Installation.txt`** in the **`slicot`** directory.

## Version v5.8

Version v5.8 of the **SLICOT** Library includes 18 new routines and 14 routines with more or few changes in the operational part of the source code, in comparison to Version v5.7 of the library.  Moreover, the (comment) lines referring to the version number have been removed in all routines, example programs, and documentation files.

**New Routines**

**`MA02AZ:`** (Conjugate) transposes all or part of a two-dimensional complex matrix.

**`MB01UY:`** Computes one of the matrix products `T := alpha*op(T)*A`, or `T := alpha*A*op(T)`, where `alpha` is a scalar, `A` is an `M-by-N` matrix, `T` is a triangular matrix, and `op(T)` is either `T` or `T'` (the transpose of `T`). A block-row/column algorithm is used, if possible. The result overwrites the array `T`.

**`MB01UZ:`** Computes one of the matrix products `T := alpha*op(T)*A`, or `T := alpha*A*op(T)`, where `alpha` is a scalar, `A` is an `M-by-N` complex matrix, `T` is a complex triangular matrix, and `op(T)` is `T`, or `T'` (the transpose of `T`), or `conj(T')` (the conjugate transpose of `T`). A block-row/column algorithm is used, if possible. The result overwrites the array `T`.

**`MB03RW:`** Solves the Sylvester equation `-A*X + X*B = C`, where `A` and `B` are complex `M-by-M` and `N-by-N` matrices, respectively, in Schur form. This routine is intended to be called only by **SLICOT** Library routine `MB03RZ`. For efficiency purposes, the computations are aborted when the absolute value of an element of `X` is greater than a given value `PMAX`.

**`MB03RZ:`** Reduces an upper triangular complex matrix `A` (Schur form) to a block-diagonal form using well-conditioned non-unitary similarity transformations. The condition numbers of the transformations used for reduction are roughly bounded by `PMAX`, where `PMAX` is a given value. The transformations are optionally postmultiplied in a given matrix `X`. The Schur form is optionally ordered, so that clustered eigenvalues are grouped in the same block.

**`MB03VW:`** Reduces the general product `A(:,:,1)^S(1) * A(:,:,2)^S(2) * ... * A(:,:,K)^S(K)` to upper Hessenberg-triangular form, where `A` is `N-by-N-by-K` and `S` is the signature array with values 1 or -1 (as exponents). The `H`-th matrix of `A` is reduced to upper Hessenberg form while the other matrices are triangularized.  Optionally, all or part of the transformation matrices are accumulated or updated.

**`SB03OS:`** Solves for `X = op(U)^H * op(U)` either the stable non-negative definite continuous-time Lyapunov equation

                    H                     2      H
               op(S) *X + X*op(S) = -scale *op(R) *op(R),

 or the convergent non-negative definite discrete-time Lyapunov equation

                    H                     2      H
               op(S) *X*op(S) - X = -scale *op(R) *op(R),

 where `op(K) = K` or `K^H` (the conjugate transpose of the matrix `K`), `S`, `R`, and `U` are `N-by-N` upper triangular matrices, and scale is an output scale factor, set less than or equal to 1 to avoid overflow in `X`.

**`SB03OZ:`** Solves for `X = op(U)^H * op(U)` either the stable non-negative definite continuous-time Lyapunov equation

                    H                     2      H
               op(A) *X + X*op(A) = -scale *op(B) *op(B),

 or the convergent non-negative definite discrete-time Lyapunov equation

                    H                     2      H
               op(A) *X*op(A) - X = -scale *op(B) *op(B),

 where `op(K) = K` or `K^H`  (the conjugate transpose of the matrix `K`), `A` is an `N-by-N` matrix, `op(B)` is an `M-by-N` matrix, `U` is an upper triangular matrix containing the Cholesky factor of the solution matrix `X`, and scale is an output scale factor, set less than or equal to 1 to avoid overflow in `X`.

**`SG03BR:`** Computes the parameters for the complex Givens rotation (`c` real, `s` complex) to annihilate the second element of a complex vector of length two. The first element of the rotated vector may be complex. This is a safer implementation of the previous **SLICOT** routine `SG03BY`. It is an adaptation for real double precision computations of the LAPACK routine `ZLARTG`.

**`SG03BS:`** Computes the Cholesky factor `U` of the matrix `X`, `X = op(U)^H * op(U)`, which is the solution of the generalized `d`-stable  discrete-time Lyapunov equation

                H            H                  2    H
               A  * X * A - E  * X * E = - SCALE  * B  * B,

 or the conjugate transposed equation

                        H            H          2        H
               A * X * A  - E * X * E  = - SCALE  * B * B ,

 respectively, where `A`, `E`, `B`, and `U` are complex `N-by-N` matrices, and `SCALE` is an output scale factor, set less than or equal to 1 to avoid overflow in `X`. The Cholesky factor `U` of the solution is computed without first finding `X`. The pencil `A - lambda * E` must be in complex generalized Schur form (`A` and `E` are upper triangular and the diagonal elements of `E` are non-negative real numbers). Moreover, it must be `d`-stable, i.e., the moduli of its eigenvalues must be less than one. `B` must be an upper triangular matrix with real non-negative entries on its main diagonal. The resulting matrix `U` is upper triangular. The entries on its main diagonal are non-negative.

**`SG03BT:`** Computes the Cholesky factor `U` of the matrix `X`, `X = op(U)^H * op(U)`, which is the solution of the generalized `c`-stable continuous-time Lyapunov equation

                H            H                  2    H
               A  * X * E + E  * X * A = - SCALE  * B  * B,

 or the conjugate transposed equation

                        H            H          2        H
               A * X * E  + E * X * A  = - SCALE  * B * B ,

 respectively, where `A`, `E`, `B`, and `U` are complex `N-by-N` matrices, and `SCALE` is an output scale factor, set less than or equal to 1 to avoid overflow in `X`. The Cholesky factor `U` of the solution is computed without first finding `X`. The pencil `A - lambda * E` must be in complex generalized Schur form (`A` and `E` are upper triangular and the diagonal elements of `E` are non-negative real numbers). Moreover, it must be `c`-stable, i.e., its eigenvalues must have negative real parts. `B` must be an upper triangular matrix with real non-negative entries on its main diagonal. The resulting matrix `U` is upper triangular. The entries on its main diagonal are non-negative.

**`SG03BZ:`** Computes the Cholesky factor `U` of the matrix `X`, `op(U)^H * op(U)`, which is the solution of either the generalized `c`-stable continuous-time Lyapunov equation

                    H                    H                      2        H
               op(A)  * X * op(E) + op(E)  * X * op(A) = - SCALE  * op(B)  * op(B),

 or the generalized `d`-stable discrete-time Lyapunov equation

                    H                    H                      2        H
               op(A)  * X * op(A) - op(E)  * X * op(E) = - SCALE  * op(B)  * op(B),

 without first finding `X` and without the need to form the matrix `op(B)^H * op(B)`. `op(K)` is either `K` or `K^H` for `K = A`, `B`, `E`, `U`. `A` and `E` are `N-by-N` matrices, `op(B)` is an `M-by-N` matrix. The resulting matrix `U` is an `N-by-N` upper triangular matrix with non-negative entries on its main diagonal. `SCALE` is an output scale factor set to avoid overflow in `U`.

**`TG01KD, TG01KZ:`** Compute for a single-input single-output descriptor system, `(A, E, B, C)`, with `E` upper triangular, a transformed system, `(Q'*A*Z, Q'*E*Z, Q'*B, C*Z)`, via an orthogonal/unitary equivalence transformation, so that `Q'*B` has only the first element nonzero and `Q'*E*Z` remains upper triangular. `TG01KZ` is the complex version.

**`TG01OA, TG01OB:`** Compute for a single-input single-output descriptor system, `(A, E, B, C)`, with `E` upper triangular, a transformed system, `(Q'*A*Z, Q'*E*Z, Q'*B, C*Z)`, via an orthogonal/unitary equivalence transformation, so that `Q'*B` has only the first element nonzero and `Q'*E*Z` remains upper triangular. `TG01OA` is the real version and `TG01OB` is the complex version. `A`, `B`, `C` are stored in an array as the block elements `(2,2)`, `(2,1)`, and `(1,2)`, respectively, and `Q` and `Z` are not accumulated. These are the main differences with **SLICOT** Library routines `TG01KD` and `TG01KZ`.

**`TG01OD, TG01OZ:`** Compute for a single-input single-output descriptor system, `(A, E, B, C)`, with `E` nonsingular, a reduced system with a "sufficiently" large feedthrough variable, using `TG01OA/TG01OB`.

**Updated Routines**

**`AB13MD:`** Computed the correct upper bound on the structured singular value for a 1-by-1 real matrix and a real uncertainty. Replaced `DFLOAT` by `DBLE` in line `834` to avoid a trouble with an Apple Silicon compiler.

**`MA02EZ:`** A new option, `SKEW = 'G'`, has been added that allows to suitably deal with the diagonal of a general square triangular matrix. This option is needed in the new routines `MB01UZ`, `SG03BS`, and `SG03BT`. Moreover, the internal loop index `J` has been modified from `2` to `I` or `I+1`, to reduce the number of cycles to the minimum values.

**`MB03RD:`** Replaced `PMAX*PMAX` by `PMAX` in the comments;  the `PMAX` value agrees to the condition number for the optimally scaled transformations used.

**`MB04BD:`** Made a correction (in comments) of the indices of the `1-by-1` or `2-by-2` quadruple diagonal blocks stored in `DWORK`. Increased `I2X2` by `1` if a `2x2` quadruple of diagonal blocks is found to have real eigenvalues (in order to check them externally). These are stored as `2x2` quadruple blocks in DWORK in that case. Made two corrections of the pointers to the locations in `DWORK` storing the quadruples with unreliable eigenvalues.

**`SB03OD:`** Many changes have been made to improve the routine, the main ones being summarized below:
  - Added code segments to control overflow. In essence, two scaling strategies are included. One strategy scales `A` and `B` if the maximum absolute value of their elements are outside a range `[SMLNUM,BIGNUM]`, where `SMLNUM = sqrt( SAFMIN )/EPS`, `BIGNUM = 1/SMLNUM`, `SAFMIN` is the safe minimum, and `EPS` is the machine accuracy. The second strategy, invoked for continuous-time equations, scales `A` and `B` if the maximum absolute values of `A` and `B` differ too much, or their minimum (maximum) is too large (small, respectively);  specifically,  this scaling is performed if `MN < MX*SMLNUM`, or `MX < SMLNUM`, or `MN > BIGNUM`, where `MN` and `MX` are the minimum and maximum, respectively, of the maximum absolute values of `A` and `B`.  Both strategies are effective and ensure the same accuracy of the results, but the second strategy reduces the number of instances when the output scaling factor, `SCALE`, is strictly smaller than `1`. Scaling of `B` is done before computing its initial `QR` or `RQ` factorization if the maximum absolute value of its elements is greater than `1/SAFMIN`; otherwise, it is done after `QR/RQ` factorization. The implementation checks out first the conditions for the second scaling strategy.
  - The auxiliary routine `SB03OU` is no longer called, and all its computations, and additional ones, are performed by `SB03OD`.  Two `QR` or `RQ` factorizations for `B` and `R*Q` or `Q'*R`, respectively (where `R` is the triangular factor of the first factorization) are done only if `M > 7*N/6`;  otherwise, a single `QR` or `RQ` factorization is used.
  - The new routine `MB01UY` is called by `SB03OD` to compute the product `R*Q` or `Q*R`, overwriting the array containing `R`;  as large block-row or block-column operations as possible (depending on the workspace length) are used.
  - More BLAS 3 `DGEMM` operations are used when updating the given `B` as `B*Q` or `Q'*B`, using again as large block-row or block-column as possible.  The previous version used BLAS 2 `DGEMV` calls if the workspace length was smaller than `N*M`.
  - The eigenvalues of `A` are computed even if `A` is given in the real Schur form, by calling the LAPACK routine `DLANV2`.  This way, the stability of `A` can be checked out also for the option `FACT = 'F'`.
  - The case when `Q` is an identity matrix is detected, to avoid its use in multiplications.
  - The postponed scaling of transformed `B` is dealt with separately for the cases `TRANS = 'N'` and `TRANS = 'T'`.
  - The minimal workspace size has been reduced by `min(N,M)`.
 The improvements allowed to reduce the computation times comparing to the previous version of this routine, sometimes by `10%`, or even `20%`.  Many changes in the comments have also been performed.

**`SB03OT:`** Deleted the lines 509-512 involving a test of the variable `TEMP` (compared to `SMIN`), which proved to be unnecessary.

**`SB03OY:`** Replaced the definition of `SMIN` in the lines 211-212 by `SMLNUM`.  The former definition set `SMIN` to a too large value for equations with elements of big magnitude, and the results were bad.

**`SG03BD:`** Many changes have been made to improve the routine, the main ones being summarized below:
  - Added code segments to control overflow. In essence, two scaling strategies are included.  One strategy scales `A`, `E`, and `B` if the maximum absolute value of their elements are outside a range `[SMLNUM,BIGNUM]`, where `SMLNUM = sqrt( SAFMIN )/EPS`, `BIGNUM = 1/SMLNUM`, `SAFMIN` is the safe minimum, and `EPS` is the machine accuracy. The second strategy scales `A`, `E`, and `B` if the maximum absolute values of `A`, `E`, and `B` differ too much, or if their minimum (maximum) is too large (small, respectively); specifically, this scaling is performed if `MN < MX*SMLNUM`, or `MX < SMLNUM`, or `MN > BIGNUM`, where `MN` and `MX` are the minimum and maximum, respectively, of the maximum absolute values of `A`, `E`, and `B`.  Both strategies are effective and ensure the same accuracy of the results, but the second strategy reduces the number of instances when the output scaling factor, `SCALE`, is strictly smaller than `1`. Scaling of `B` is done before computing its initial `QR` or `RQ` factorization if the maximum absolute value of its elements is greater than `1/SAFMIN`; otherwise, it is done after `QR/RQ` factorization.  The implementation checks out first the conditions for the second scaling strategy.  The scaling factors of `E` may be set equal to those for `A`, to preserve stability in the discrete-time case.
  - Modified the sequences for the transformation of the right hand side, to use BLAS 3 `DGEMM` calls (possibly in a loop, if workspace is not large enough).
  - Used the new routine, `MB01UY`, for transforming the solution back by block-row or block-column operations. The result overwrites the array storing the upper triangular part of the reduced equation solution.
 - The call of `DGGES` in the sequence to find the optimal workspace is omitted if `A` and `E` are given in generalized Schur form.
 - Avoided the back transformation operations with the matrices `Q` and `Z`, if they are identity.
 - Made some additional comments on the output contents of `A`, `E`, `Q`, and `Z` arrays, as well as on the input contents of `Q` and `Z` arrays when `FACT = 'F'`.

**`SG03BU:`** Safer computation of the 1-by-1 diagonal blocks of the solution. Replaced loops of `DSCAL` calls by calls to the LAPACK routine `DLASCL`. Replaced the calls to BLAS 1 routine `DROTG` by calls to the safer LAPACK routine `DLARTG`.

**`SG03BV:`** Safer computation of the 1-by-1 diagonal blocks of the solution, to avoid overflows for badly scaled equations. Replaced loops of `DSCAL` calls by calls to the LAPACK routine `DLASCL`. Replaced the calls to BLAS 1 routine `DROTG` by calls to the safer LAPACK routine `DLARTG`.

**`SG03BW:`** Replaced the calls to the **SLICOT** routines `MB02UV` and `MB02UU` by calls to the equivalent LAPACK routines `DGETC2` and `DGESC2`, respectively. Replaced loops of `DSCAL` calls by calls to the LAPACK routine `DLASCL`. Replaced some loops of assignment statements by calls to `DCOPY`.

**`SG03BX:`** The code has practically been rewritten. Although the source code of the new version is somewhat longer than of the previous version, the object code was reduced to about a half of the initial one.
  - Replaced all `2-by-2` local arrays by the needed scalars, and replaced the calls to the BLAS routines `DGEMM` and `DGEMV` involving these arrays by the necessary scalar operations, made in a compact manner. In particular, all rotations are applied directly, updating only the needed elements. When possible, two rotations are applied at once (from the left and right). Also, only the necessary elements of the matrices M1 and M2 of the reduced equations are now computed.
  - Replaced the calls to the **SLICOT** routine `SG03BY` by calls to the new, safer routine `SG03BR`.
  - Added a sequence of code to standardize the matrix `E`, i.e., to become diagonal with `E(1,1)` non-negative; then recomputed the shift.
  - Safer computation of the `1-by-1` diagonal blocks of the solution, to avoid overflows for badly scaled equations.
  - Used simpler formulas and computations when numerically appropriate, but more sophisticated ones are invoked when needed.
  - When `E(1,1)` is not very small for discrete-time equation, an alternative approach is employed for computing the trailing diagonal element of the solution; this approach factors a rank-one Hermitian matrix by solving a special symmetric eigenvalue problem. The eigenvector corresponding to eigenvalue `1` is quickly found.

**`TB01MD, TB01ND:`** Moved the sequence for initializing `U` before the Quick Return section, to set `U` to identity also for the case `M = 0` or `P = 0`, respectively.

**`TB01ND:`** Added the condition `P <= N`, since the observer Hessenberg form is not defined for `P > N`.

## Version v5.7

This is the first version uploaded to GitHub under a BSD-3-Clause license.