File Release4.History ===================== This file contains a short description of the previous modifications performed in the Release 4.0 of SLICOT Library since November 1999, in reverse chronological order. The latest changes are listed in the file Release.Notes. Both files will be updated from time to time. ======================================================================== Update Jan. 05, 2005 Corrected Bugs: ============== AB13DD : Used SAFMAX for TM if the pencil has too large eigenvalues. Moved the statements for "Return if the lower bound is zero" just before "Start the modified gamma iteration for the Bruinsma-Steinbuch algorithm". Set INFO = 4 if ITER > MAXIT. AB13DX : Set the minimal workspace IWORK and CWORK to 0 and 1 if P = 0 or M = 0. IB03AD : Replaced INFO by INFOL in the two calls of MD03AD in the loop 10. SB10ED : Replaced A by A - I in the call of SB10PD, since SB10PD performs the corresponding tests for continuous-time systems. Array A is modified internally, but it is restored on exit. Made the Quick return conditions more restrictive, as needed. In some cases, slightly reduced the size of IWORK. Improvements: ============ SB08ED, : Replaced the scalar argument ALPHA by an array argument of SB08FD length 2. The first entry contains the stability degree, and the second entry contains the stability margin. Previously, these entities have been equal. Documentation: ============= AB13DD : LCWORK can be 1 also if P = 0 or M = 0. AB13DX : Added that A and B are changed (as previously mentioned) only if M > 0 and P > 0. If P = 0 or M = 0, the lengths of IWORK and CWORK may be 0 or 1, respectively. SB08ED, : Updated for the changes in SB08ED and SB08FD (see above). SB08FD SB10ED : Replaced A by A - I in two formulas, and mentioned that the array A is modified internally, but it is restored on exit. In some cases, slightly reduced the size of IWORK. Updated test programs: ===================== TSB08ED, : Updated for the changes in SB08ED and SB08FD. TSB08FD Updated MEX-files: ================= aresolc : set SEP = 0 if N = 0. deadbeat : arrays V and IWORK are also allocated (minimally) for ISTAIR = -1 and for ISTAIR = -2 and ISTAIR = 1, respectively; array U is initialized to identity matrix if WithU = 0 and ISTAIR = 0 or 1; set INFO = 0 if SB06ND is not called (i.e., ISTAIR < 0); printed INFO1, if INFO1 <> 0. fstoep : set INFO = 0 if UPDATE is TRUE, but MB02DD is not called (i.e., M = 0). gsystra : set M = 0 and P = 0 for task = 3 and 2, respectively. Deallocated IWORK if task = 5. mucomp : replaced maximally allocated arrays with arrays allocated exactly as needed. The former code, containing arrays with sizes specified in PARAMETER statements did not work under Linux with g95 Fortran compiler and MATLAB 14. Removed the arrays DK and DT and used X instead. Mentioned the limitations concerning the size of a real block and the sum of sizes. syscf : updated for the changes in SB08ED and SB08FD. Updated M-files: =============== cmpoles : replaced "cmppoles" (in comments) by "cmpoles". gsystra : replaced "decresingly" (in comments) by "decreasingly". mucomp : mentioned the limitations concerning the size of a real block and the sum of sizes. syscf : updated for the changes in SB08ED and SB08FD. test_dlsdz : used tol = tole/10 for Unix machines, instead of 0. Updated demonstrations M-files: ============================== basicdemo : inserted a final command uiwait to force the user to close this demo before exploring another demo. Simplified the callback to the button "close". cmp_red : changed the current directory to SLdemos (in order to use the data files stored there), and finally restored the initial working directory. Saved the variable "saved" also when vars is empty. Contents : replaced "slicot_demos" by "slicotdemos". dex_wident, : replaced the statement dex_widentc if exist( 'range' ) ~= 1 | isempty( range ), by the code sequence if exist( 'range' ) ~= 1, range = 40; elseif isempty( range ), since "range" is already used in the Stats toolbox. slicot_demos : the file was deleted, but its contents have been incorporated in slicotdemos (moved from subdirectory SLTools in SLdemos). slicotdemos : added the functionality of slicot_demos. Added two tuning parameters (pause_duration and reproducible) so that the user can adapt more easily the demos. slmodred : deleted the reference to cmp_red, since now cmp_red is called from slredemo. slredemo : added a new demo function, comparing SLICOT and MATLAB model reduction functions. slwidemo5 : set n = size(sys.a,1); after calling slmoen4. Replaced the file SEED.mat by Seed.mat. ======================================================================== Update April 14, 2004 Note: All updates have been saved in the compressed library files, slicot.tar.gz and slicotPC.zip. Corrected Bugs: ============== MB01RD, : In the Quick return section, R is set also if N = 0. MB01RU, MB01RX MB02CU : RNK is set to 0 when K = 0 and TYPEG = 'D'. Array A2 is referenced instead of A1 at the second DNRM2 call in loop 50. MB02CV : Replaced F1(LEN,J) by F1(J,LEN) in line 639. MB02FD : Added the term 4*K in the formula for the workspace length. MB02ID : Corrected the formula of optimal workspace in line 486 and performed slight code improvements. MB02JD : Added the parameter tests P >= 0 and S >= 0. MB02JX : Defined DWORK(2:3) and RNK for quick returns. MB02KD : Before the Quick return section, C is set to 0, if BETA = 0. Quick return is also performed if R = 0. TB01PD : In the Quick return section, IWORK is set to 0. TG01AD : Replaced B(I,J) by B(I,1) in statement 379. TG01JD : Moved the initialization of INFRED before the test for Quick return, to be used in all cases. Several calls to DCOPY (most often with negative increments) have been replaced with in-line FOR loops, to avoid possible overwritting when using optimized BLAS libraries. Specifically, such changes have been performed in the following routines and in the specified loops: AB05MD, : Loops 10, 20, 30, and 40. AB05PD, AB05QD, AB05RD AB05ND : Loops 10, 30, and 40. FB01SD : Loop 20. FB01TD : Loop 30. IB01MD : Loops 10, 380, and 390. IB01ND : Loops 20 and 30. Specifically, Z is set to an identity matrix if JOBZ = 'I', and Z and TAU are set to zero matrix and vector, respectively, if JOBF = 'I'. IB03AD, : DCOPY before the first call of TF01MX. IB03BD MB02DD : Loop 40. MB02MD : The call of DCOPY. MB02ND : Loops 360 and 380. MB05MD, : The call of DCOPY. MB05MY MD03BD : The last call of DCOPY. NF01BP : The last call of DCOPY (in the line 641). NF01BS : Loops 50, 100 (second DCOPY call) and 120 (second DCOPY call). SB01BD : All two calls of DCOPY. SB01BX : All three calls of DCOPY. SB03OT : The call of DCOPY in the line 477. SG03BD : The call of DCOPY in the line 691 is made only if N > 1. TB01VD : Loop 40. The call of DCOPY in the line 355 is made only if L <> 2. TB01VY : Loop 10. TB04BV : Loop 70. TB04BW : Loop 50. All Fortran and Matlab tests have been redone using a library version built with the debug option activated. Several routines, listed below, have been changed in order to avoid exceeding array bounds. While this could previously happen, no writting action was actually performed, but just wrong references to array elements in the calls. The changes also enabled to make slight improvements of some codes (AB05ND, FB01RD, IB01MD, IB01MY, MB01UW, SB02MU), see below. AB01ND : The last call of DLASET (in line 425) is made only if N >= IQR. AB05MD : Quick return is now made if MAX( N, MIN( M1, P2 ) ) = 0. The calls of DLACPY and DLASET containing references to subarrays are made if the corresponding dimensions are nonzero. AB05OD : Quick return is now made if MAX( N, MIN( M, P1 ) ) = 0. The calls of DLACPY and DGEMM containing references to subarrays are made if the corresponding dimensions are nonzero. AB05PD, : Quick return is now made if MAX( MAX( N, MIN( M, P ) ) = 0. AB05QD The calls of DLACPY and DLASET containing references to subarrays are made if the corresponding dimensions are nonzero. AB07MD : Quick return is now made if MAX( N, MINMP ) = 0. The loops 10 and 20 are entered only if N > 0, and the loop 30 is entered only if JOBD = 'D' and MIN( M, P) > 0. AB08ND : The loop labelled 240 is rerun only if DINFZ > 0. IB01RD : The call of DCOPY in line 439 is made only if N > 1. MA02CD : N-1 was replaced by N-2 in the expressions for the counter I in the loops 10 and 20. MB02CD : The calls of MB02CY in lines 340, 397, 488, and 546 are done only if N > I. The calls of DTRSM in lines 277 and 425, and the calls of DLACPY in lines 287, 380, 435, and 529 are made only if N > 1. MB02CU : All calls of DLARFG, DLARF, DROT, DAXPY, and DSCAL are now made such that array bounds are not exceeded. MB02CV : Defined LEN = MAX( N - K, 0 ), instead of LEN = N - K. MB02CX : All calls of DLARFG and DLARF are now made such that array bounds are not exceeded. MB02CY : The two calls of DLASET are made only if Q > 1. MB02DD : The first call of MB02CY in each of the loops 30 and 60 is made only if M > I. Also, the calls to DLACPY before the loops 10 and 40 are made only if N > 1. MB02ED : The first call of DGEMM in each of the loops 10 and 20 is made only if N > I. Also, when needed, certain calls to DGEMM, DLACPY, and DTRSM are made only if N > 1. MB02FD : The two calls of DTRSM are made only if N > 1. MB02GD : The two calls of DTRSM are made only if NL > 1. The call of DLASET in the line 340 is made only if I > 1. MB02HD : The call of DLASET in the line 364 is made only if LENR < N*L and I > 1. MB02JX : The calls of DLASET, DLACPY, MA02AD, and MB02KD are now made such that array bounds are not exceeded. MB02ND : The call of DLASET in the line 504 is made only if NL > 1. MB02UV : Before the loop 100, DSCAL and DGER are called only if N > 1. MB03VY : The call of DORGQR in line 183 is made only if NH > 1. TB01XD : The transposition of D was moved before the Quick return section, so that D is modified even if N = 0. N-1 was replaced by N-2 in the expressions for the counter J in the loops 20 and 30. TB01YD : The call of DSWAP after the loop 10 is made if N is odd and N > 2. Note: In many cases, there were no indication of errors, but the mentioned changes increase the reliability. Slight improvements: =================== AB05ND : The loop 20 (containing a call of DCOPY with negative increments) has been replaced with a call to DLACPY from a copy of the matrix C. The IF THEN ELSE clauses have been modified accordingly. FB01RD : The negative increments in the DCOPY call in loop 20 have been replaced by 1. IB01MD : The counter J of loop 10 now starts with 2. IB01MY : The counter I of loop 10 now starts with 2, and the negative increments in the DCOPY call have been replaced by 1. MB01UW : The cases M = 1 and SIDE = 'L' or N = 1 and SIDE = 'R' are now solved by calling DSCAL and RETURN. SB02MU : The loops 280 and 320 have been rewritten, by separating the cases UPLO = 'U' and UPLO = 'L'; each case has now its own loop, executed N-1 times, and the remaining DCOPY call is done outside the loop. Documentation: ============= MB02CV, : Minor changes. MB02CX, MB02GD, MB02HD, MB02JD, MB02JX, MB02KD MB02FD : Added the term 4*K in the formula for the workspace length. TG01ID : Replaced the first occurrence of E(1,1) in (1) by E(1,2), and similarly replaced the first occurrence of A(1,1) in (3) by A(1,2). TG01JD : Replaced Hessenebrg by Hessenberg in line 90. Updated mex-files: ================= gsyscom.f, : The calls of DCOPY with negative increments have been invert.f, replaced by suitable FOR loops. sysconn.f, TotalLS.f Updated m-files: =============== Contents.m : added descriptions for the new mexfiles and m-files. test_persch : modified the tests when only the eigenvalues between ilo and ihi are computed (for given argument index). test_polezero : modified the tests for defining Znrm in lines 772, 773, and those for Pnrm in lines 835, 836, to avoid the warning "Future versions will return empty for empty == scalar comparisons" in MATLAB version 6.1. New mexfiles: ============ muHopt.f : computes the mu optimal or H_inf controller. New m-files: =========== muHopt.m : help function for the mexfile muHopt.f. slihinf.m : computes the H_inf optimal controller given a state space model. slimju.m : computes the mu optimal controller given a state space model. Also, computes the mu norm of the closed loop system. test_slimju.m : test for slimju.m. ======================================================================== Update January 3, 2004 Note: All updates have been saved in the compressed library files slicotPC.zip and slicot.tar.gz. Moreover, the files update.zip and update.tar.gz contain the updated and new files. Corrected Bugs: ============== AB01ND : Set Z, TAU, and DWORK(1) in the Quick return part if N > 0. Specifically, Z is set to an identity matrix if JOBZ = 'I', and Z and TAU are set to zero matrix and vector, respectively, if JOBF = 'I'. AB01OD : Set U, V, and DWORK(1) in the Quick return part, if needed. AB08MD : Set DINFZ to 0 before Quick return, and made a quick return if MIN( M, P ) is 0, in addition to N = 0. As consequences, the calls of DLACPY and TB01ID are made only if NN and PP are nonzero. The loop labelled 240 is entered only if N > 0. AB08ND : Replaced DWORK(NP*M+NP+1) by DWORK(NP*M+N+1) in the call of TB01ID. AB09JV, : The calls of DTGSYL and DTRSYL are now made only if N > 0. AB09JW FB01QD : A Quick return is now made if N = 0. FB01RD : A Quick return is now made if N = 0. The routine now also works for P > N. Replaced N-P-1 by MAX( N-P-1, 0 ) in the MB04JD call. FB01SD, : A Quick return is now made if MAX( N, P ) = 0, and a return FB01TD is made before applying the transformations to the last column of the pre-array, if N = 0. X is also suitably set when JOBX = 'N'. FB01TD : The routine now also works for M > N. Replaced N-MP1 by MAX( N-MP1, 0 ) in the MB04ID call. FB01VD : A Quick return is now made if MAX( N, L ) = 0. MB02MD : Quick return is now done if min( M, N+L ) = 0, and, if M = 0, C is set to an identity matrix and X to a zero matrix. Another return is added after the loop 40 if L = 0. Also, the WHILE loop 80 is entered only if RANK < min( M, N+L ). MB02ND : Quick return is now done if min( M, N+L ) = 0, and, if M = 0, C is set to an identity matrix, X to a zero matrix, and INUL is set to .TRUE.. The definition of the variable SUFWRK was slightly changed. Another return is made before starting STEP 4 if L = 0. MB03MD : Added the statement IF ( N.EQ.1 ) RETURN after the line calling the function MB03MY. SB08CD, : Replaced the arguments involving NQ-1 in the call of TB01XD by SB08ED MAX( 0, NQ-1 ). SB10MD : Replaced LORD by ORD in the second call of TB05AD (19-th argument). SB10YD : Replaced TOL by TOLL in the two calls of DGELSY. SG03BX : Replaced the in-line calculation of the eigenvalues of the 2-by-2 pencil A - lambda * E by a call to LAPACK Library routine DLAG2. Fixed a bug in the normalization of the (2,2) entry of the upper triangular factor in the QR factorization of B. TB01UD : Set FNRM to ONE if FNRM < TOLDEF. Z and TAU are now set when M = 0, but N > 0 (Quick return). Replaced MAX( 1, ITAU-1 ) by ITAU-1 in the call of DORGQR, for getting a correct transformation matrix when B is negligible. TB01ZD : Set TAU(1) = ZERO if N = 0 (line 324). TB05AD : The Quick return is now done if N = 0 and G is then set to the complex zero matrix when MIN( M, P ) > 0. DSWAP and DSCAL are now called only if M > 0 or P > 0. TG01HX : Set SVLMAX to ONE if SVLMAX < RCOND. Improvements: ============ FB01RD : The routine now also works for P > N. FB01TD : The routine now also works for M > N. MB02ND : The necessary workspace length was somewhat reduced. SB06ND : The code was slightly changed to allow the computation of the minimum norm feedback matrix also when the uncontrollable part of the system is zero. Documentation: ============= AB08ND : Replaced N+1 by MAX(N,M)+1 and MAX(N,P)+1 for the dimensions of the arrays KRONR and KRONL, respectively, to allow dealing with with systems with more inputs/outputs than states. FB01QD, : Updated the statement about the content of K when JOBK = 'N', FB01RD or JOBK = 'K' and INFO = 1. FB01SD, : Added the term "innovation" in the statement about the FB01TD content of QINV on exit. FB01VD : Replaced "discrete algebraic Riccati" by "discrete-time difference Riccati". MB02MD : The columns N+L-RANK : N+L of C on exit contain the transformed right singular vectors of C, modified as described in Step 3 of the TLS algorithm, if RANK > 0 and L > 0. Otherwise, the right singular vectors are not transformed. MB02ND : The right singular vectors corresponding to the small singular values of C are transformed as described in Step 4 of the PTLS and L > 0. Otherwise, the right singular vectors are not transformed. The formula for the necessary workspace length was somewhat reduced. Updated m-files: =============== Contents.m : added descriptions for the new mexfiles and m-files. New mexfiles: ============ cfsys.f : constructs the state-space representation of a system from the factors of its left or right coprime factorization. deadbeat.f : constructs the minimum norm feedback matrix performing "deadbeat control" on an (A,B)-pair. Kfiltupd.f : computes a combined measurement and time update of one iteration of the Kalman filter. polezero.f : computes the normal rank, poles, zeros, and the Kronecker structure of the system pencil for a standard or descriptor system. specfact.f : computes the spectral factorization of a real polynomial. TotalLS.f : solves the Total Least Squares (TLS) problem using a singular value decomposition (SVD) approach or a Partial SVD (PSVD) approach. New m-files: =========== cfsys.m : help function for the mexfile cfsys.f. cf2ss.m : constructs the state-space representation of a system from the factors of its left or right coprime factorization. test_cfsys.m : test for cfsys.f and cf2ss.m. deadbeat.m : help function for the mexfile deadbeat.f. ckstair.m : checks that a system is in a staircase form. deadbt.m : constructs the minimum norm feedback matrix performing "deadbeat control" on a system. test_ckstair.m : test for ckstair.m. test_deadbeat.m: test for deadbeat.f and deadbt.m. Kfiltupd.m : help function for the mexfile Kfiltupd.f. tLSfilt.m : computes a combined measurement and time update of the recursive least-squares filter. It is intended for testing the gateway Kfiltupd for task = 3. test_Kfiltupd.m: test for Kfiltupd.f. polezero.m : help function for the mexfile polezero.f. nrank.m : computes the normal rank of the transfer-function matrix of a standard system. polzer.m : computes the normal rank, poles, zeros, and the Kronecker structure of the system pencil for a standard or descriptor system. slpole.m : computes the poles of a standard or descriptor system. slzero.m : computes the normal rank, zeros, and the Kronecker structure of the system pencil for a standard or descriptor system. test_polezero.m: test for polezero.f and the associated M-files. specfact.m : help function for the mexfile specfact.f. polysf.m : computes the spectral factorization of a real polynomial. test_specfact.m: test for specfact.f and polysf.m. TotalLS.m : help function for the mexfile deadbeat.f TLS.m : solves the TLS problem using a SVD or partial SVD approach. housh.m : computes a Householder transformation which modifies one element of a vector and annihilates all other elements. test_TotalLS.m : test for TotalLS.f and associated M-files. ======================================================================== Update September 5, 2003 Note: All updates have been saved in the compressed library files, slicot.tar.gz and slicotPC.zip. Corrected Bugs: ============== AB05MD, : Added new tests for dealing with matrices with zero dimensions, AB05ND, and to ensure that LDC (or LDCx) is 1 when N (or Nx) is 0. AB05OD, Details are given below. AB05PD AB05MD : Replaced C(I1,1) by C(1,I1) in a DGEMM call (in the line 491). Replaced I1 = N2 + 1 by I1 = MIN( N2 + 1, N ), and similarly for I2. Added tests IF ( N1.GT.0 ) when forming the off-diagonal blocks of matrix A, and when forming C. AB05ND : For OVER = 'O', increased the workspace by N1*P1 (to save C1, which otherwise is overwritten by part of C), and replaced M1*M1 by M1*(M1+1), if M1 > N*N2. The BLAS 2 code in the loop labelled 30 is used for M1 <= N*N2, and another loop with additional M1 workspace is used for M1 > N*N2. Replaced the second ONE by ZERO in the DGEMV call before the label 30. The statements containing references to rows/columns indexed N1+1 are not executed if N2 = 0. Added tests IF ( N1.GT.0 ), IF ( N2.GT.0 ), IF ( M1.GT.0 ), or IF ( P1.GT.0 ) whenever needed. Added some code to avoid overwritting parts of A and B if OVER = 'O'. AB05OD, Added the test N.GT.0 for the DLASCL call for C. AB05PD AB05RD : The two DGEMM calls for computing Ac and C1 using K are executed only for BETA.NE.0 and N.GT.0. The DGEMM call for computing Cc is executed only for N.GT.0. AB05SD : The first DGETRS call is executed only if N.GT.0. After the calculation of Cc and Dc, the routine returns if N = 0. The last DGEMM call is executed only if JOBD = 'D'. AB13DD : The bugs in applying the permutations for the balancing transformations on B and C have been fixed. In particular, LAPACK routines DGGBAK and DGEBAK cannot be used. MB05MD, : The backward transformation of the eigenvectors for MB05MY balancing (when BALANC = 'S') was moved from MB05MY to MB05MD, since it must be done after premultiplication with the orthogonal matrix performing the reduction to real Schur form. The scaling is also used when computing the inverse of the matrix V. The scaling factors are stored in DWORK(2), ..., DWORK(N+1). MB05OD : RERL is set to 0 when SIZE <= EMNORM. SB01BD : Added few suitable statements after the call of SB01BY, when a simple 1x1 block is uncontrollable. SB01BY : If R is 0, it is reset to EPSM, the machine accuracy, to avoid a division by 0. If Y in the loop 10 is 0, then the loop is exited, for the same reason. SB08CD, : Replaced MIN( N, M, P ).EQ.0 by MIN( N, P ).EQ.0 in the SB08ED Quick return part. That is, if M = 0, but MIN( N, P ) > 0, the calculation continues appropriately, delivering the result for a system (A,C). SB08DD, : Replaced MIN( N, M, P ).EQ.0 by MIN( N, M ).EQ.0 in the SB08FD Quick return part. That is, if P = 0, but MIN( N, M ) > 0, the calculation continues appropriately, delivering the result for a system (A,B). SB08FD : INFO on exit from SB01BY is now checked, since, due to roundoff, an uncontrollable 2x2 block could appear, which actually corresponds to a double real eigenvalue; one eigenvalue is then eliminated, and the procedure restarted. Such a block wasn't previously removed. SB10LD : replaced NP2 by M2 in the penultimate call of DGEMM (bug found by A. Markovski). SG03BX : Added the statement INFO = 0. (Proposed by Klaus Schnepper.) TF01MD, : The computation of the output vectors when n = 0 (non-dynamic TF01ND system case) was made correct. TG01AD : LSCALE or RSCALE is set to 1 if L > 0 but N = 0, or L = 0 but N > 0, respectively. (Quick return sequence.) TG01CD : The call of DORMQR for updating B is made only if M > 0. TG01AD, : The calling statements containing indexed references to the TG01ED, arrays B and/or C are now not executed if M = 0 and/or P = 0, TG01FD, so that B and/or C can have 0 length. TG01HX TG01ID : Added the statement IF ( P.EQ.0 .OR. NR.EQ.0 ) LBE = MAX( 0, N - 1 ) before the last call of TB01XD. TG01JD : The leading dimensions N, M, and P in the calls of DLACPY have been replaced by MAX(1,N), MAX(1,M), and MAX(1,P), respectively. The leading dimension LDC has been replaced by MAX( LDC, P ) in the call of TG01AD and the first two calls of TG01HX, and by MAX( LDC, M ) in the last two calls of TG01HX. This is needed for N = 0. INFRED has been set for Quick return. Improvements: ============ MB03RD : - Two new options, SORT = 'C' and SORT = 'B' have been added, which allow a "closest-neighbour" strategy to be used for selecting the block to be added to the current block for block diagonalization. The previous options SORT = 'N' and SORT = 'S' use the "closest to the mean" strategy. The new strategy often performs block-diagonalization when the old strategy does not succeed, and/or produces better conditioned transformations. The closest-neighbour strategy was suggested by Pascal Gahinet. - Quick return (for N = 0) has been included. New Routines: ============ SB10AD : to compute the matrices of an H-infinity optimal n-state controller, using modified Glover's and Doyle's 1988 formulas. SB10MD : to perform the D-step in the D-K iteration (continuous-time case). SB10YD : to fit a supplied frequency response data with a stable, minimum phase SISO (single-input single-output) system represented by its matrices A, B, C, D. SB10ZP : to transform a SISO system [A,B;C,D] by mirroring its unstable poles and zeros in the boundary of the stability domain, thus preserving the frequency response of the system, but making it stable and minimum phase. Documentation: ============= AB05ND : For OVER = 'O', increased the workspace by N1*P1, and replaced M1*M1 by M1*(M1+1), if M1 > N*N2. MB03RD : Two new options, SORT = 'C' and SORT = 'B' have been added, which allow a "closest-neighbour" strategy to be used for selecting the block to be added to the current block. MB05MD, : Few changes mainly due to moving the backward transformation MB05MY of the eigenvectors for balancing (when BALANC = 'S'), see above. TB05AD : Replaced "If BALEIG = 'B' or 'E' or BALEIG = 'A'" by "INITA = 'G'" in the description of A, B, and C (on exit). TG01JD : Replaced "INDRED(7)" by "INFRED(7)" in the description of IWORK. The title texts for Phases 1 to 4 (METHOD) now mention that "all" corresponding eigenvalues are eliminated. libindex, : added links to the new document files. support Updated M-files: =============== Contents.m : added descriptions for the new MEX-files and M-files. Mentioned that sparse matrix format is not currently supported. slH2norm.m, : replaced the test "if Ts > 0" by "if Ts ~= 0", to slHknorm.m, allow dealing with discrete-time ss systems with slinorm.m unspecified sampling time. test_aresol.m, : deleted the commands "clear all" since it appears to test_aresolc.m, suspend the execution of MATLAB 6.1, and to avoid test_conred.m, destroying the user data. test_genleq.m, test_linmeq.m, test_syscf.m, test_syscom.m, test_sysred.m, test_sysred_tools.m, test_systra.m test_datana.m, : deleted the commands "clear variables" to avoid test_Hessol.m, destroying the user data. test_Hnorm.m test_dsim.m : set "warning off" before the loop for comparison with lsim, to suppress warnings produced by MATLAB 6.5 \toolbox\control\ctrlguis\rfinputs and lsim, and set "warning on" after the loop. test_linorm.m : replaced the command "E = eye(n);" by "E = [];" for testing this new option (for standard systems). Updated MEX-files: ================= linorm.f : Standard systems (with E an identity matrix) may now be defined by setting E as an empty matrix (i.e., with 0 rows and/or columns), to reduce the memory requirements. syscf.f : The routine names SB08AD and SB08BD in two error messages have been replaced with the correct ones (SB08ED and SB08FD, respectively). New MEX-files: ============ bldiag.f : performs block-diagonalization of a general matrix or a matrix in real Schur form. condis.f : performs a transformation on the parameters (A,B,C,D) of a system, which is equivalent to a bilinear transformation of the corresponding transfer function matrix. gsyscom.f : transforms a descriptor system, by equivalence transformations, to a controllable or observable staircase form, or to a reduced (controllable, observable, or irreducible) form. gsystra.f : performs various equivalence transformations for descriptor systems with scaling, generalized Schur form, etc. invert.f : computes the dual or inverse of a linear (descriptor) system. ldsimt.f : computes the output response of a linear discrete-time system. The input and output trajectories are stored column-wise (each column contains all inputs or outputs measured at a certain time instant). slmexp.f : computes the matrix exponential and optionally its integral. sysconn.f : computes a state-space model (A,B,C,D) for various inter-connections of two systems given in state-space form. sysfconn.f : computes, for a given state-space system (A,B,C,D), the closed-loop system (Ac,Bc,Cc,Dc) corresponding to the output, or mixed output and state, feedback control law. New M-files: =========== bldiag.m : help function for the MEX-file bldiag.f. bdiag.m : block diagonalization of a general matrix or a matrix in real Schur form. test_bldiag.m : test for bldiag.f and bdiag.m. condis.m : help function for the MEX-file condis.f. slc2d.m : bilinear transformation of a continuous-time system to a discrete-time system. sld2c.m : bilinear transformation of a discrete-time system to a continuous-time system. test_condis.m : test for condis.f and associated M-files. gsyscom.m : help function for the MEX-file gsyscom.f. slgconf.m : reduces a descriptor system to controllable staircase form. slgobsf.m : reduces a descriptor system to observable staircase form. slgminr.m : reduces a descriptor system to an irreducible form. test_gsyscom.m : test for gsyscom.f and associated M-files. gsystra.m : help function for the MEX-file gsystra.f. slgsbal.m : balances the system matrix for a descriptor system. slgsHes.m : transforms the pair (A,E) of a descriptor system to a generalized Hessenberg form. slgsQRQ.m : transforms the pair (A,E) of a descriptor system to a QR- or RQ-coordinate form. slgsrsf.m : transforms the pair (A,E) of a descriptor system to a real generalized Schur form. slgsSVD.m : transforms the pair (A,E) of a descriptor system to a singular value decomposition (SVD) or SVD-like coordinate form. test_gsystra.m : test for gsystra.f and associated M-files. invert.m : help function for the MEX-file invert.f. sldual.m : computes the dual of a standard system. slinv.m : computes the inverse of a standard system. test_invert.m : test for invert.f and associated M-files. ldsimt.m : help function for the MEX-file ldsimt.f. dsimt.m : computes the output response of a linear discrete-time system. The trajectories are stored column-wise. test_ldsimt.m : test for ldsimt.f and dsimt.m. slmexp.m : help function for the MEX-file slmexp.f. slexpe.m : computes the matrix exponential basically using an eigenvalue/eigenvector decomposition technique, but a diagonal Pade approximant with scaling and squaring, if the matrix appears to be defective. slexpm.m : computes the matrix exponential using a diagonal Pade approximant with scaling and squaring. slexpi.m : computes the matrix exponential and optionally its integral, using a Pade approximation of the integral. test_slmexp.m : test for slmexp.f and associated M-files. sysconn.m : help function for the MEX-file sysconn.f. slapp.m : appends two systems in state-space form. slfeed.m : feedback inter-connection of two systems in state-space form. slpar.m : parallel inter-connection of two systems in state-space form. slser.m : series inter-connection of two systems in state-space form. slspar.m : rowwise inter-connection of two systems in state-space form. test_sysconn.m : test for sysconn.f and associated M-files. sysfconn.m : help function for the MEX-file sysfconn.f. slosfeed.m : closed-loop system for a mixed output and state feedback control law. slofeed.m : closed-loop system for an output feedback control law. test_sysfconn.m: test for sysfconn.f and associated M-files. ======================================================================== Update March 27, 2003 Note: All updates have been saved in the compressed library files, slicot.tar.gz and slicotPC.zip. Corrected Bugs: ============== AB13BD : recomputed the Frobenius norm of D if NR.GT.0. AB13ED : set LOW and HIGH also for N = 0. AB13FD : set BETA and OMEGA also for N = 0. DF01MD : works for N >= 5 (not 3, as assumed before). MA02HD : replaced IF ( J.NE.M ) THEN by IF ( J.LT.M ) THEN before the loop labelled 60. MB03SD : removed the values JOBSCL = 'P' and 'B', for which the Hessenberg form of the matrix A'' would not be preserved. Annihilated the lower triangle below the Hessenberg form when JOBSCL = 'S'. Improved efficiency by using BLAS 3 calls to DSYMM as much as possible. SB03QD, : replaced the first two occurrences of IF( UPDATE ) THEN by SB03SD IF( NOFACT .OR. UPDATE ) THEN. SB10SD : redefined IWRK = IWV + N*N before the second call of SB02SD. Improvements: ============ SB02MD : the scaling factor used is now returned in DWORK(2). (The reciprocal condition number of the matrix A is returned in DWORK(3) in the discrete-time case.) SB02OD : - Internal scaling has been incorporated, to increase accuracy and reliability for poorly scaled equations. - The natural tendency of the QZ algorithm to get the largest eigenvalues in the leading part of the matrix pencil is exploited for discrete-time Riccati equations, by computing the unstable eigenvalues of the permuted matrix pencil. - A standard eigenproblem is solved for continuous-time equations when G is given (JOBB = 'G'). SB02OY : the second matrix is not built (and memory is saved) for continuous-time problems with G given and identity E matrix, since a standard eigenproblem is solved in this case. SB02RD : if JOB = 'X', the scaling factor used is now returned in SEP. SG02AD : the natural tendency of the QZ algorithm to get the largest eigenvalues in the leading part of the matrix pencil is exploited for discrete-time Riccati equations, by computing the unstable eigenvalues of the permuted matrix pencil. Documentation: ============= DF01MD : works for N >= 5 (not 3, as assumed before). MB03SD : removed the values JOBSCL = 'P' and 'B', for which the Hessenberg Schur form of the matrix A' would not be preserved. SB02MD : the scaling factor used is now returned in DWORK(2). (The reciprocal condition number of the matrix A is returned in DWORK(3) in the discrete-time case.) SB02OD : if JOBB = 'B', mentioned that the arrays Q, R, and L are modified internally, but are restored on exit. DWORK(3) returns the scaling factor used internally. SB02OY : the second matrix is not built (and memory is saved) for continuous-time problems with G given and identity E matrix, since a standard eigenproblem is solved in this case. SB02RD : the scaling factor used is now returned in SEP. SG02AD : DWORK(4) is always set to the scaling factor (or to 1, when scaling is not used). DAESolver,: added documentation files for these SLICOT interfaces to ODESolver, various solvers. FSQP, KINSOL libindex : added links to the new document files. Updated m-files: =============== Contents.m : added descriptions for the new mexfiles and m-files. slcaregs.m, : updated X2 for the scaling factor used. slcares.m, slcaresc.m, sldaregs.m, sldaregsv.m, sldares.m, sldaresc.m test_aresol.m, : updated some tests for the scaling factor used. Also, non-zero test_aresolc.m coupling matrix L is now used. Updated mex-files: ================= aresol.f, : updated for changes in SB02MD, SB02OD, and/or SB02RD. aresolc.f The scaling factor used is now returned when a basis for the stable subspace is computed. Updated files for MATLAB 6.5: ============================ bst.m, : files on the left have been updated by changing bta.m, a command like btabal.m, discr = sys.ts > 0; btabal_cf.m, to bta_cf.m, discr = double(sys.ts > 0); cfconred.m, since MATLAB 6.5 assumes discr to be a logical array, fwbconred.m, and this is not accepted by the called mex-files. fwbred.m, These changes are also valid for MATLAB 5.3 and 6.1. fwhna.m, hna.m, lcf.m, lcfid.m, rcf.m, rcfid.m, spa.m, spabal.m, spabal_cf.m, spa_cf.m, test_conred.m All SLICOT test functions for SLICOT mex-files and m-files have been run under MATLAB 6.5, and the changed funtions re-checked under MATLAB 5.3 and 6.1. Executable mex-files for MATLAB 6.5 on Windows platforms have been generated and made available on the SLICOT ftp site (file dllfiles65.zip, subdirectory SLTools). These mex-files are significantly smaller than the corresponding MATLAB 6.1 mex-files. Updated zip-files: ================= basic_mex.zip : updated for the changes performed. modred_mex.zip, data_mex.zip : new .zip file for data analysis calculations. conred_mex.zip : now includes help functions for bstred, conred, fwehna, fwered, and sfored mex-files. Replaced the files bstred.dll and conred.mat from the archive in subdirectory SLToolboxes5 by their Matlab 5.3 versions. Now, the included tests are also working under Matlab 5.3. hinf_mex.zip : updated for the changes performed (updated dishin, dlsdp, dlsdz). New mexfiles: ============ arecond.f : estimates condition and forward errors for Lyapunov and algebraic Riccati equations. garesol.f : solves descriptor algebraic Riccati equations. datana.f : performs various transforms of real or complex vectors (data analysis calculations). Hessol.f : performs analysis and solution of a system of linear equations with an upper Hessenberg coefficient matrix. Hnorm.f : computes various system norms (Hankel norm, H2 norm) and complex stability radius. Hamileig.f : computes the eigenvalues of a Hamiltonian matrix using the square-reduced approach. New m-files: =========== arecond.m : help function for the the mexfile arecond.f. carecond.m : reciprocal condition number of a continuous-time algebraic Riccati equation. darecond.m : reciprocal condition number of a discrete-time algebraic Riccati equation. lyapcond.m : reciprocal condition number of a Lyapunov equation. steicond.m : reciprocal condition number of a Stein equation. test_arecond.m : test for arecond.f and associated M-files. garesol.m : help function for the the mexfile garesol.f. slgcare.m : solution of a descriptor continuous-time algebraic Riccati equation. slgdare.m : solution of a descriptor discrete-time algebraic Riccati equation. test_garesol.m : test for garesol.f and associated M-files. cndclyap.m, : auxiliary test files for test_garesol.m. cnddlyap.m, cndricc.m, cndricd.m datana.m : help function for the the mexfile datana.f. sincos.m : sine or cosine transform of a real signal. slDFT.m : discrete Fourier transform of a signal. sliDFT.m : inverse discrete Fourier transform of a signal. slHart.m : discrete Hartley transform of a real signal. slconv.m : convolution of two real signals using either FFT or Hartley transform. sldeconv.m : deconvolution of two real signals using either FFT or Hartley transform. slwindow.m : anti-aliasing window applied to a real signal. test_datana.m : test for datana.f (via the M-files). sncs.m : auxiliary test file for sincos.m. Hessol.m : help function for the the mexfile Hessol.f. Hesscond.m : computes an estimate of the reciprocal of the condition number of an upper Hessenberg matrix. Hessl.m : solves a set of systems of linear equations with an upper Hessenberg coefficient matrix. test_Hessol.m : test for Hessl.m and Hessol.f. Hnorm.m : help function for the the mexfile Hnorm.f. slH2norm.m : computes H2/L2 norm of a system. slHknorm.m : computes Hankel-norm of a stable projection of a system. slstabr.m : computes complex stability radius. test_Hnorm.m : test for slH2norm.m, slHknorm.m, slstabr.m, and Hnorm.f. Hamileig.m : help function for the the mexfile Hamileig.f. Hameig.m : computes the eigenvalues of a Hamiltonian matrix. test_Hameig.m : test for Hameig.m and Hamileig.f. bstred.m, : help functions for the corresponding mexfiles. conred.m, fwehna.m, fwered.m, sfored.m New contributed files: ===================== munorm.m : mu-norm computation file wrote by Craig T. Lawrence and Andre L. Tits. test_scalar.m, : test files for munorm. test_mixed.m, test_complex.m robpole.m, : robust pole assignment codes by Y. Yang and complex_pair.m, Andre L. Tits. The underlying algorithm is somewhat one_column.m, faster than the algorithm due to Kautsky, Nichols and real_pair.m Van Dooren (1985), and it can handle complex prescribed eigenvalues. ======================================================================== Update October 1, 2002 Corrected Bugs: ============== MB02CU : the calls to DLARF and DROT in the loops 100, 120, 130, and 150 are made only if K.GT.J. Similarly, the calls to DSCAL and DAXPY in the loops 130 and 150. MB02HD : the conditions on ML and NU were set in a way so that the first MIN(M*K,N*L) columns of T generically have full rank. Checked the added conditions on ML and NU. Replaced "MIN( I-1, N )*L," by "MIN( I-1, N-1 )*L," in the first DGEMM call in the loop 20. Replaced "LENR.LT.N*L" by "LENR.LT.MIN( ML+NU+1,N )*L" in the test of the loop 40. Added " .AND. J.GT.1" in the test of the loop 80. MB02ID : replaced N*L by M*K (twice) in the DO 10 statement. Replaced "LDB" by "LDC" in the second DGEMM call after the label 40. The calls to DGEMM for "Block Gaussian elimination to B" and "Block Gaussian elimination to C" are made only if LEN.GT.L. Replaced "NB" by "RB" and "RC", respectively, in the two calls of DTRMM. MB02JD : The three calls to MB02KD are made only if N.GT.1. The first MB02CV call is made only if LEN.GT.KK, and the next DLACPY call is made only if M.GT.1. Replaced the check of the bound on LDR (see section "Documentation"). MB02JX : replaced "5*L" by "9*L" (twice) in the value of LDWORK. The first two calls to MB02KD are made only if N.GT.1. The DLACPY call after the label 70 is made only if N.GT.2. Replaced "(M-1)*K" by "M*K" in a DLASET call after 90. The first MB02CV call is made only if LEN.GT.KK. The fifth argument in the two MB02CV calls was changed to PP. The DLACPY call in the loop 170 is made only if LEN.GT.2*L. Documentation: ============= MB02HD : deleted the condition "M*K >= L" in the description of M. Added new conditions in the description of ML and NU. MB02JD : replaced MAX by MIN in the description of parameter R on exit, the previous formula being MAX( N*L,(N-P+1)*L )-by-MIN( S*L, MIN( M*K,N*L )-P*L ) and replaced accordingly the bound on LDR. MB02JX : replaced "5*L" by "9*L" (twice) in the value of LDWORK. SB10ZD : added the document file for this routine. SG02AD : added comments related to the solution of the associated optimal control problem when matrix R is singular. libindex : added links to the new document files. On-line HTML documentation files for SLICOT routines performing model reduction on parallel computers, using either direct or iterative methods have been made available on the SLICOT ftp site. Updated m-files: =============== Contents.m : added descriptions for the new mexfiles and m-files. pass.m : if DISCR = 1, replaced ALPHA = Inf with ALPHA = 0 (in the help part). persch.m : replaced % [W,Z] = PERSCH(A,JOB,COMPZ,INDEX) with JOB = 0, or by % [AO,Z] = PERSCH(A,JOB,COMPZ,INDEX) with JOB = 0, or (in the help part). polass.m : if discr = 1, replaced alpha = Inf with alpha = 0 (second item in section Comments). Updated mex-files: ================= perschur.f : annihilated the elements of A(:,:,2:p) below the main diagonal and of A(:,:,1) below the first subdiagonal also for JOB = 0. polass.f : if discr = 1, replaced alpha = Inf with alpha = 0 (second item in section Comments). Also, replaced the use of DWORK for storing the number of unmodified, assigned, and uncontrollable poles by TMP (needed for N = M = 0). New mexfiles: ============ fstoeq.f : computes an orthogonal-triangular/trapezoidal decomposition (with partial column pivoting) of a (banded) block Toeplitz matrix. Executable Matlab 6.1 mexfiles have been generated for Sun Sparc-Solaris and are stored in the compressed file mexsolfiles6.tar.gz, subdirectory MatlabTools/Unix/SLTools of the SLICOT ftp site. The mexfiles fstoeq, perschur and polass are not included. The available mexfiles have been tested on a Sun sparc SUNW Ultra-2 (machine hardware sun4u, OS version 5.8). New m-files: =========== SB10ZD.m : script called by test_dlsdz.m for testing SB10ZD. btoeplitz2.m : generates general block Toeplitz matrices. fstlsq.m : solves linear least-squares problems min(B-T*X) or finds the minimum norm solution of T'*Y = C where T is a (block) Toeplitz matrix with full column rank, given the first block column and the first block row of T. fstmul.m : computes the matrix-vector products X = T*B for a block Toeplitz matrix T, given the first block column and the first block row of T. fstqr.m : computes the orthogonal-triangular decomposition of a (block) Toeplitz matrix T and solves associated linear least-squares problems, given the first block column and the first block row of T, and assuming that the first min(size(T)) columns of T have full rank. test_fstoeq.m : test script for the fast solver for linear least-squares problems involving an mxk-by-nxl (banded) block Toeplitz matrix. test_pass.m : script for testing the m-file pass and mexfile polass. test_persch.m : script for testing the m-file persch and mexfile perschur. New demo files: ============== wident_demo.zip : contains all needed files to demonstrate some features of the toolbox for Wiener system identification. All previously available Matlab 5.3 demo files have been moved in a subdirectory called SLdemos5, and new versions for Matlab 6.1 replaced the existing files in the subdirectory SLdemos. New directory: ============= prebuilt : contains prebuilt SLICOT object library files. The following files are currently included: slicot_solaris77.tar.gz : containing slicot.a - object library file for Sun sparc Ultra-2 machines, built using Sun WorkShop Compiler FORTRAN 77 5.0. slicot_solaris90.tar.gz : containing slicot.a - object library file for Sun sparc Ultra-2 machines, built using Sun WorkShop Compiler FORTRAN 90 2.0. slicot_win32.zip : containing slicot.lib - object library file for PC Windows 9x/NT/ME/2000 machines, built using Compaq Visual Fortran 6.5. slicot_linux77.tar.gz : containing slicot.a - object library file for Red Hat Linux 7.2 machines, built using Fortran 77. ======================================================================== Update June 29, 2002 Note: All updates have been saved in the compressed library files, slicot.tar.gz and slicotPC.zip. Corrected Bugs: ============== AB13FD : replaced the statement PARAMETER ( CONE = ( ONE, ZERO ), $ RTMONE = ( ZERO, ONE ) ) by PARAMETER ( CONE = ( 1.0D0, 0.0D0 ), $ RTMONE = ( 0.0D0, 1.0D0 ) ) in order to avoid compilation errors on certain Fortran compilers. IB03AD, : replaced the statement IB03BD IWARN = MAX( IWARN, 10*IWARNL ) by IF ( IWARN.GT.100 ) THEN IWARN = MAX( IWARN, ( IWARN/100 )*100 + 10*IWARNL ) ELSE IWARN = MAX( IWARN, 10*IWARNL ) END IF after calling the Levenberg-Marquardt routine(s) for initialization of the nonlinear part. MB03NY : replaced the statement PARAMETER ( CONE = ( ONE, ZERO ) ) by PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) in order to avoid compilation errors on certain Fortran compilers. MC01OD : moved the statements REP(1) = ONE IMP(1) = ZERO (previously placed before the DO statement) before the return for K = 0. So, REP(1) and IMP(1) are now initialized also when K = 0. This is not a bug, but more convenient in usage. MC01PD : moved the statement P(1) = ONE (previously placed before the commented WHILE statement) before the return for K = 0. So, P(1) is now initialized also when K = 0. This is not a bug, but more convenient in usage. SB02OD : replaced the statements WRKOPT = U(1,1) IF ( LJOBB ) RCONDL = U(2,1) by WRKOPT = DWORK(1) IF ( LJOBB ) RCONDL = DWORK(2) after the call of SB02OY. New Routines: ============ MC01PY : to compute the coefficients of a real polynomial P(x) from its zeros. The coefficients are stored in decreasing order of the powers of x. SG02AD : to solve either the continuous-time or discrete-time algebraic Riccati equations for descriptor systems. TB04BD : to compute the transfer function matrix G of a state-space representation (A,B,C,D) of a linear time-invariant multivariable system, using the pole-zeros method. Each element of the transfer function matrix is returned in a cancelled, minimal form, with numerator and denominator polynomials stored either in increasing or decreasing order of the powers of the indeterminate. TB04BV : to separate the strictly proper part from the constant part of a proper transfer function matrix. TB04BW : to compute the sum of an P-by-M rational matrix and a real P-by-M matrix. TB04BX : to compute the gain of a single-input single-output linear system, given its state-space representation (A,b,c,d), and its poles and zeros. The matrix A is assumed to be in an upper Hessenberg form. TB04CD : to compute the transfer function matrix G of a state-space representation (A,B,C,D) of a linear time-invariant multivariable system, using the pole-zeros method. The transfer function matrix is returned in a minimal pole-zero-gain form. New test programs: ================= TSG02AD : test program for SG02AD.f. TTB04BD : test program for TB04BD.f. TTB04CD : test program for TB04CD.f. Documentation: ============= MC01OD : replaced the phrase "If K = 0, these arrays are not referenced." (in the description of the arguments REP, IMP) by "If K = 0, then REP(1) is set to one and IMP(1) is set to zero.". MC01PD : replaced the phrase "If K = 0, this array is not referenced." (in the description of the argument P) by "If K = 0, then P(1) is set to one.". SB02OD : replaced "output" by "state" in the description of the argument Q. Deleted the obsolete paragraph referring to LDWORK in section FURTHER COMMENTS. * : for user's and programmer's callable newly added routines. libindex, : added links to the new document files. support support added links to MB02WD.html and MB02XD.html. Updated m-files: =============== Contents.m : included the description of the new mexfiles and m-files (mainly for the Wiener System Identification Toolbox). test_genleq : made the tolerance for slgesg 10^-8 (in the Basic Toolbox for standard and generalized state space systems and transfer matrix factorizations). Updated NICONET web page files: ============================== The files describing Tasks I.A, III.A, III.B, and IV.A have been updated. New mexfiles: ============ polass.f : performs (partial) pole assignment. perschur.f : computes the periodic Hessenberg or periodic Schur decomposition of a matrix product. New m-files: =========== pass.m : performs (partial) pole assignment (calling polass). persch.m : computes the periodic Hessenberg or periodic Schur decomposition of a matrix product (calling perschur). ======================================================================== Update April 9, 2002 Corrected Bugs: ============== AG08BD : - replaced the test ELSE IF( TOL.GT.ONE ) THEN by ELSE IF( TOL.GE.ONE ) THEN (also in the documentation). - changed the dimension of the integer arrays KRONR, INFE, and KRONL to (N+M+1), (1+MIN(L+P,N+M)), and (L+P+1), respectively. AG08BY : inserted the statement PR = 0 in the "Quick return" section. IB03BD : - replaced the statement CALL TF01MX( N, M, L, NSMP, DWORK(AC), LDAC, U, LDU, DWORK(IX), $ DWORK(Z), NSMP, DWORK(JWORK), LDWORK-JWORK+1, $ INFO ) by CALL DCOPY( N, DWORK(IX), 1, X(NTHS+1), 1 ) CALL TF01MX( N, M, L, NSMP, DWORK(AC), LDAC, U, LDU, X(NTHS+1), $ DWORK(Z), NSMP, DWORK(JWORK), LDWORK-JWORK+1, $ INFO ) - replaced the statement WORK(2) = WORK(2) + DWORK(JWORK+1) by WORK(2) = MAX( WORK(2), DWORK(JWORK+1) ) after the call of MD03BD for the initialization step. - corrected the value of JWORK (and hence of LDWORK) for the whole optimization. - corrected the returned value of IWARN. MB04VX : replaced throughout the variable CE by CJE, and replaced the variable CA by CJA just before the loop 80, and inside that loop. MD03BF : replaced INF0 by INFO. SB03MX : the two calls of DSYMV in the loop 120 are now made only if L2 < N. SB03OD : the strict lower triangle of the resulted Cholesky factor U is set to zero. Also, if M = 0, then U is set to zero. These are not bugs, but simplify the external work, and ensure compatibility with SG03BD. In addition, LDWORK may be 1 for M = 0. SG03BD : the call of DLASET before the "DO 140" statement is now made only if N > 1. TB01UD : the test "IF( IQR.LE.N )" is made before calling DLASET for annihilating the trailing blocks of B. TB01VD, : replaced the statements TB01VY IF ( LAPPLY ) (THEN) by IF ( LAPPLY .AND. TI.NE.ZERO ) (THEN) TB01VY : replaced the statements IF ( N.EQ.0 ) $ RETURN by IF ( N.EQ.0 ) THEN RETURN ELSE IF ( L.EQ.0 ) THEN CALL DCOPY( N, THETA(N*M+1), 1, X0, 1 ) RETURN END IF in the Quick Return part of the code. Corrected Bugs in Example Programs: ================================== TAG08BD : changed the dimension of the integer arrays KRONR, INFE, and KRONL to (N+M+1), (1+MIN(L+P,N+M)), and (L+P+1), respectively. Corrected Bugs in m/mexfiles: ============================ plot_ye.m, : inserted commas before the calls of "num2str" in "ylabel" plot_yu.m statements. (These files are contained in the SLIDENT demonstration file slident_demo.zip, which was updated.) ex_wident.m: simplified the logic when initx = 1; initialized errL, errN, errLm, and errNm to 0. test_wiener: corrected the value of ldwin. Updated Documentation: ===================== AG08BD : changed the declared dimension of the integer arrays KRONR, INFE, and KRONL to (N+M+1), (1+MIN(L+P,N+M)), and (L+P+1), respectively. IB03BD : - replaced ML by L + M in the definition of LTHS (for LDWORK). - the possible values for i and j, defining INFO, can be from 1 to 7 (separated the values for FNC with IFLAG = 1 and 2). NF01AY : inserted the missing section heading "ARGUMENTS". SB03OD : the strict lower triangle of the resulted Cholesky factor is set to zero. LDWORK may be 1 for M = 0. libindex, : added links to the new document files. support New Routines: ============ IB03AD : to compute a set of parameters for approximating a Wiener system in a least-squares sense, using a neural network approach and a conjugate gradients-based Levenberg-Marquardt algorithm. MB02WD : to solve a system of linear equations Ax = b, with A symmetric, positive definite, or, in the implicit form, f(A, x) = b, where y = f(A, x) is a symmetric positive definite linear mapping from x to y, using the conjugate gradient algorithm without preconditioning. MB02XD : to solve a set of systems of linear equations, A'*A*X = B, or, in the implicit form, f(A)*X = B, with A'*A positive definite, using symmetric Gaussian elimination. MD03AD : to find the parameters PARA for a function F(X, PARA) that give the best approximation for Y = F(X, PARA) in a least-squares sense. NF01BA : template routine to evaluate the functions and Jacobian matrices for optimizing the parameters of the nonlinear part of a Wiener system (initialization phase) using SLICOT Library routines MD03AD and IB03AD. NF01BB : template routine to evaluate the functions and Jacobian matrices for solving a standard nonlinear least squares problem using SLICOT Library routines MD03AD and IB03AD. NF01BU : to compute the matrix J'*J + c*I, for the Jacobian J given in a compressed form. NF01BV : to compute the matrix J'*J + c*I, for the Jacobian J fully given, for one output variable. NF01BW : to compute the matrix-vector product x <-- (J'*J + c*I)*x, where J is given in a compressed form. NF01BX : to compute x <-- (A'*A + c*I)*x, where A is an m-by-n real matrix, c is a scalar. New test programs: ================= TIB03AD : test program for IB03AD.f. TIB03BD : test program for IB03BD.f. TMD03AD : test program for MD03AD.f. New directories: =============== plicmr : contains the index of the available parallel SLICOT library routines (for large order model reduction). plicmr/doc : contains the documentation html files for the available parallel SLICOT library routines (for large order model reduction). Documentation: ============= * : for all user's and programmer's callable newly added routines listed above. New mexfiles: ============ widentc.f : computes a discrete-time model of a Wiener system using a neural network approach and a Levenberg- Marquardt algorithm with a Cholesky-based or conjugate gradients solver. New m-files: =========== ex_widentc.m : script demonstrating the use of the mexfile widentc. Updated Toolbox: =============== Wident : to incorporate the updates and the new routines mentioned above for identifying Wiener systems. ======================================================================== Update February 22, 2002 Note: All updates (including those in October 3, 2001) have been saved in the compressed library files, slicot.tar.gz and slicotPC.zip. Corrected Bugs: ============== Most of the changes below have been performed to initialize some variables in certain cases. Some of them are related to the optimal workspace length. AB01MD,: added "B1 = B(1)" when N = 1. TB01ZD AB09HD : replaced IWARNL by IWARN in the call of AB09IX. AB09HX : permuted the values for the workspace pointers KU and KV, so that the space for V (with pointer KV) could be reused. KTAU is now initialized to KV (not to 1) before the block of code containing the loop 70, to avoid possible overwriting of the workspace and/or possible non-initialization. AB09KD : replaced (twice) DWORK(1) = MAX( WRKOPT, DWORK(1) ) by WRKOPT = MAX( WRKOPT, DWORK(1) ) and added DWORK(1) = MAX( WRKOPT, DWORK(1) ) before return. BB01AD : replaced the statement "Q(ISYMM) = 100.0D0" by "Q(ISYMM) = 10.0D0", for example 3.1. BB02AD : replaced the statement "IF ((NR(2) .EQ. 3) .OR. (NR(2) .EQ. 6)) THEN" by "IF ((NR(2) .EQ. 5) .OR. (NR(2) .EQ. 6)) THEN" (for example 1.5), and the statement "C(3,11) = ONE" by "C(4,11) = ONE" (for example 1.11). Deleted the second statement setting TEMP for example 2.4. Also, set X for example 2.1. MA02CD : the three calls to DSWAP are made only if I1.GT.0. MB01PD : moved the statement "ISUM = 0" before "IF( NBL.GT.0 ) THEN", and the tests of the argument TYPE before those of the argument SCUN, to avoid possibly non-initialized variables ISUM and ITYPE. MB03PY : replaced "IF ( RANK.LT.M ) THEN" by "IF ( RANK.LT.K ) THEN". MB03WD : moved the statements HP00 = ONE HP01 = ZERO before "IF( K.GT.L+1 ) THEN", in the loop labelled 80, to avoid possibly non-initialized variables HP00 and HP01. SB02RD : inserted the statement IF ( .NOT.JOBA ) $ WRKOPT = 0 after IF ( .NOT.JOBX ) THEN to initialize the variable WRKOPT for JOB = 'C' or JOB = 'E'. SB03OD : - inserted ELSE WRKOPT = 0 at the end of the "IF ( NOFACT ) THEN" block. - put the update of WRKOPT just after the first calls of DGERQF and DGEQRF. SB10RD : inserted LWAMAX = 0 before the "IF( ND1.GT.0 ) THEN" block. TB01LD : inserted WRKOPT = 0 in the "ELSE" part of the "IF( LSAME( JOBA, 'G' ) ) THEN" block. TG01ED : inserted IR1 = 1 in the body of the "IF( LN2.EQ.0 ) THEN" block. Improved Routines: ================= AB09HY : replaced BLAS Level 2 routines DGERQ2 and DORGR2 with corresponding BLAS Level 3 routines. Corrected Bugs in Example Programs: ================================== TAB09MD, TAB09ND : in the computation of LDWORK, replaced the lines $ MAX( NMAX, MMAX, PMAX) + 5 + $ ( NMAX*(NMAX+1)/2) ) ) by $ MAX( NMAX, MMAX, PMAX ) + 5 ) + $ ( NMAX*(NMAX+1) )/2 ) TAB13MD : added the lines * .. Intrinsic Functions .. INTRINSIC MAX TBB01AD, In format statements, replaced "rows of B" by "columns of B" TBB02AD : and "columns of C" by "rows of C". These changes also resulted in BB01AD.res and BB02AD.res. TBB02AD : replaced "N" by "P" in the loop 60 (twice). About 50 compressed tar files "name".tar.gz, calling the above routines or example programs, have been updated on the ftp site. Corrected Bugs in benchmark data files: ====================================== BB012091.dat, : inserted a space between two consecutive data values. BB012092.dat, BB02113.dat BB012091.dat, : deleted a ; character from the end of a line. BB012092.dat Corrected Bugs in m/mexfiles: ============================ syscom.f : defined the default tolerance for single-input (m = 1) or single-output (p = 1) systems, and decided nonminimality for such systems if the Frobenius norm of B (if m = 1), or the Frobenius norm of C (if p = 1), is smaller than the given or default tolerance. findBD.f : the line below, which exceeded 72 characters, WRITE( TEXT, '(''A must have '',I6,'' rows and columns'')' ) N was splitted in two lines. New Routines: ============ IB03BD : to compute a set of parameters for approximating a Wiener system in a least-squares sense, using a neural network approach and a MINPACK-like Levenberg-Marquardt algorithm. MB02YD : to solve a system of linear equations A*x = b , D*x = 0 , in the least squares sense, with D is a diagonal matrix, given a QR factorization with column pivoting of A. MB04OW : to perform the QR factorization ( U ) = Q*( R ), where U = ( U1 U2 ), R = ( R1 R2 ), ( x' ) ( 0 ) ( 0 T ) ( 0 R3 ) where U and R are (m+n)-by-(m+n) upper triangular matrices, x is an m+n element vector, U1 is m-by-m, T is n-by-n, stored separately, and Q is an (m+n+1)-by-(m+n+1) orthogonal matrix. MD03BA : interface to SLICOT Library routine MD03BX, for solving standard nonlinear least squares problems using SLICOT routine MD03BD. MD03BB : interface to SLICOT Library routine MD03BY, for solving standard nonlinear least squares problems using SLICOT routine MD03BD. MD03BD : to find the parameters PARA for a function F(X, PARA) that give the best approximation for Y = F(X, PARA) in a least-squares sense, using a MINPACK-like Levenberg-Marquardt algorithm. MD03BF : template routine to evaluate the functions and Jacobian matrices for solving a standard nonlinear least squares problem using SLICOT Library routine MD03BD. MD03BX : to compute the QR factorization with column pivoting of an m-by-n matrix J (m >= n), that is, J*P = Q*R, where Q is a matrix with orthogonal columns, P a permutation matrix, and R an upper trapezoidal matrix with diagonal elements of nonincreasing magnitude, and to apply the transformation Q' on the error vector e (in-situ). The 1-norm of the scaled gradient is also returned. Standard nonlinear least squares problems are solved. MD03BY : to determine a value for the parameter PAR such that if x solves the system A*x = b , sqrt(PAR)*D*x = 0 , in the least squares sense, where A is an m-by-n matrix, D is an n-by-n nonsingular diagonal matrix, and b is an m-vector, and if DELTA is a positive number, DXNORM is the Euclidean norm of D*x, then either PAR is zero and ( DXNORM - DELTA ) .LE. 0.1*DELTA , or PAR is positive and ABS( DXNORM - DELTA ) .LE. 0.1*DELTA . It is assumed that a QR factorization, with column pivoting, of A is available. Standard nonlinear least squares problems are solved. NF01AD : to compute the output of a Wiener system. NF01AY : to compute the output of a set of neural networks. NF01BD : to compute the Jacobian of a Wiener system. NF01BE : template routine to evaluate the functions and Jacobian matrices for optimizing the parameters of the nonlinear part of a Wiener system (initialization phase), using SLICOT Library routine MD03BD. NF01BF : template routine to evaluate the functions and Jacobian matrices for optimizing all parameters of a Wiener system, using SLICOT Library routine MD03BD. NF01BP : to determine a value for the Levenberg-Marquardt parameter PAR such that if x solves the system J*x = b , sqrt(PAR)*D*x = 0 , in the least squares sense, where J is an m-by-n matrix, D is an n-by-n nonsingular diagonal matrix, and b is an m-vector, and if DELTA is a positive number, DXNORM is the Euclidean norm of D*x, then either PAR is zero and ( DXNORM - DELTA ) .LE. 0.1*DELTA , or PAR is positive and ABS( DXNORM - DELTA ) .LE. 0.1*DELTA . The matrix J is the current Jacobian matrix of a nonlinear least squares problem, provided in a compressed form by SLICOT Library routine NF01BD. It is assumed that a block QR factorization, with column pivoting, of J is available. NF01BQ : to determine a vector x which solves the system of linear equations J*x = b , D*x = 0 , in the least squares sense, where J is an m-by-n matrix, D is an n-by-n diagonal matrix, and b is an m-vector. The matrix J is the current Jacobian of a nonlinear least squares problem, provided in a compressed form by SLICOT Library routine NF01BD. It is assumed that a block QR factorization, with column pivoting, of J is available. NF01BR : to solve one of the systems of linear equations R*x = b , or R'*x = b , in the least squares sense, where R is an n-by-n block upper triangular matrix, with a block diagonal structure, plus an additional block column to the right, with the upper triangular diagonal submatrices R_k square, and of the same order, except for the last. The diagonal elements of each block R_k have nonincreasing magnitude. The matrix R is stored in a compressed form, as returned by SLICOT Library routine NF01BS. If the matrix R has not a full rank, then a least squares solution is obtained. Optionally, the transpose of the matrix R can be stored. NF01BS : to compute the QR factorization of the Jacobian matrix J, as received in compressed form from SLICOT Library routine NF01BD, and to apply the transformation Q on the error vector e (in-situ). The factorization is J*P = Q*R, where Q is a matrix with orthogonal columns, P a permutation matrix, and R an upper trapezoidal matrix with diagonal elements of nonincreasing magnitude for each block column. The 1-norm of the scaled gradient is also returned. NF01BY : to compute the Jacobian of the error function for a neural network (for one output variable). TB01VD : to convert the linear discrete-time system given as (A, B, C, D), with initial state x0, into the output normal form, with parameter vector THETA. The matrix A is assumed to be stable. The matrices A, B, C, D and the vector x0 are transformed, so that on exit they correspond to the system defined by THETA. TB01VY : to convert the linear discrete-time system given as its output normal form, with parameter vector THETA, into the state-space representation (A, B, C, D), with the initial state x0. TF01MX : to compute the output sequence of a linear time-invariant open-loop system given by its discrete-time state-space model with an (N+P)-by-(N+M) general system matrix S (the input and output trajectories are stored differently from SLICOT Library routine TF01MD). TF01MY : to compute the output sequence of a linear time-invariant open-loop system given by its discrete-time state-space model (A,B,C,D), where A is an N-by-N general matrix (the input and output trajectories are stored differently from SLICOT Library routine TF01MD). New test programs: ================= TMD03BD : test program for MD03BD.f. Test programs for F77: ===================== All test (example) programs which contained MAX and/or MIN intrinsic functions calls in PARAMETER statements have now a version without these calls, in order to be compliant with the Fortran 77 standard. The codes have not been optimized: the arguments of those intrinsic functions have just been added together, and the function references removed. The modified files (over 100) are stored in the subdirectory examples77. For convenience, this subdirectory contains all the example programs, data (.dat) and results (.res) files. New directory: ============= The Matlab 5.3 toolboxes have been saved in a new subdirectory, called "SLToolboxes5", of the MatlabTools directory of the ftp site. The former subdirectory "SLToolboxes" now contains the Matlab 6 versions of all files, including .dll files. The Matlab 5.3 files will not be updated in the future. Documentation: ============= AB09HX : deleted "*HNORM(A,B,C)" and the text referring to it in the description of NR. added " and ORDSEL = 'A'" in the last line of the description of TOL2. changed the numbering of the formulas in the section METHOD. AB09HY : replaced "Aw = A - Bw*Cw" by "Aw = A - Bw*Hw" in the section PURPOSE, and added the mention that matrix A must be stable and in a real Schur form. AB09KD : added M in the two MAX formulas for LIWORK, and made other minor changes. AB13AD, : added "(if INFO = 0 )" in the definition of the function AB13AX, value. AB13BD, MB03NY AB13CD : added "If INFO = 0," in the definition of the function value. AB13MD : replaced "cummulative" by "cumulative" (in the internal comments). MA02GD : replaced " If IPIV is negative" by " If INCX is negative" in the definition of INCX. * : for programmer's callable newly added routines. Updated m-files: ================ fstdemo5.m : replaced "cummulative" by "cumulative". New mexfiles: ============= arebench.f : generates ARE examples. New m-files: ============ aredata.m : generates ARE examples. New toolbox: =========== Wident : beta-version of the toolbox for nonlinear, Wiener system identification will be posted by the end of February 2002. ======================================================================== Update October 3, 2001 Note: ==== The changes described below are not yet performed in the compressed library files (slicot.tar.gz and slicotPC.zip), but they are included in the files update.tar.gz (for Unix platforms) and update.zip (for Windows platforms). Corrected Bugs: ============== AB09IY : replaced the test after the label 10 by IF( ALPHAO.LT.ONE .AND. NV.GT.0 ) THEN and the test after the label 50, and the next line, by IF( ALPHAC.LT.ONE .AND. NW.GT.0 ) THEN KTAU = N*NNW + 1 Also, the matrix D is not referenced if JOB = 'B' or 'F', hence LDD may be 1 in this case. IB01PD : replaced "TOLL" by "TOLL1" in the second test involving RCOND1, and in the second call of MB03OD. IB01BD,: replaced "ELSE IF( LDWORK.GE.1 ) THEN" by "ELSE". IB01PD, IB01PX, IB01PY SB16AD : the workspace formula includes N, when WEIGHT = 'N' and EQUIL = 'S'. SB16BD : the workspace formula for ORDSEL = 'F' and NCR = N has been added in the code. Also, some minor changes have been performed in the comments (corrected "facorization", etc.). Corrected Bugs in example programs: ================================== TAB09ID : updated the formula for LDWORK. Corrected Bugs in m/mexfiles: ============================= slinorm.m : replaced the test "elseif size(fpeak0) == 1" by the safer test "elseif isinf(fpeak0)". aresol.f, : put the line "IP = 5" in front of the statement "IF ( NRHS.GE.6 ) THEN". aresolc.f findBD.f : inserted "LDWMIN = 2" in the IF branch setting "LDWORK = 2". find_models.m (in ident_demo.zip) : replaced "rcnds = zeros(k,14);" by "rcnds = zeros(k,12);" in the first branch "if nout == 3", and "rcnds(k,:) = rcnd;" by "rcnds(k,:) = rcnd';", in the second such if branch. New Routines: ============ SB10ZD to compute the matrices of a positive feedback controller in the Discrete-Time Loop Shaping Design Procedure. New test programs: ================= TSB10ZD : test program for SB10ZD.f. Documentation Updates: ===================== AB09ID : inserted a line "LRCF = 0, and" before "LMINR = 0, ..." if WEIGHT = 'L' or 'N' or NW = 0, in the formulas for LDWORK. Other minor changes have also been performed. AB09IX : replaced S by R in the description of the parameter SCALEO. SB16AD : the workspace formula includes N, when WEIGHT = 'N' and EQUIL = 'S'; the reference to non-existing bibliographic item 4 has been replaced to item 2. SB16CD : deleted one of the two consecutive "of" in the description of JOBCF. New mexfiles: ============= dlsdz.f : Loop Shaping Design of discrete-time systems. New m-files: ============ Contents.m : lists the SLICOT m-files and mexfiles and their function. dlsdz.m : help function for the Loop Shaping Design routine for discrete-time systems. test_dlsdz.m associated test function for dlsdz.m. Parallel Codes: ============== A new file, plicmr.tar.gz, has been added to the SLICOT root directory. It contains 3 model reduction routines, 7 associated routines (solvers for coupled stable Lyapunov/Stein equations, stable Sylvester equations, etc.), and 12 lower-level routines, all for parallel computers. Details are given in the incorporated .html documentation, accessible from the file index.html, included in the archive. ======================================================================== Update June 29, 2001 Removed Bugs: ============ AB08NX : the call of DLASET setting ABCD(NU+SIGMA+1,SIGMA) is made only if RO1 > 1. AB09CX : replaced the statement "NKR1 = NR1 + KR" by the correct one "NKR1 = MIN( NMINR, NR1 + KR )". MB04ZD : a symplectic rotation is now applied on rows/columns J+2 of A and Q/G only for J < N-1. SB16AY : Replaced the ELSE statement before the comment line "Define B1 and C2 for Ge21" by "ELSE IF( PERF ) THEN". New Routines: ============ AB08MD to compute the normal rank of the transfer-function matrix of a state-space model (A,B,C,D). AB09JD : to compute a reduced order model (Ar,Br,Cr,Dr) for an original state-space representation (A,B,C,D) by using the frequency weighted optimal Hankel-norm approximation method. AB09JV : to construct a state-space representation (A,Bs,Cs,Ds) of the projection of V*G or conj(V)*G containing the poles of G, from the state-space representations (A,B,C,D) and (AV-lambda*EV,BV,CV,DV), of the transfer-function matrices G and V, respectively. G is assumed to be a stable transfer-function matrix and the state matrix A must be in a real Schur form. AB09JW : to construct a state-space representation (A,Bs,Cs,Ds) of the projection of G*W or G*conj(W) containing the poles of G, from the state-space representations (A,B,C,D) and (AW-lambda*EW,BW,CW,DW), of the transfer-function matrices G and W, respectively. G is assumed to be a stable transfer-function matrix and the state matrix A must be in a real Schur form. AB09JX to check stability/antistability of finite eigenvalues with respect to a given stability domain. AB13DD to compute the L-infinity norm of a continuous-time or discrete-time system, either standard or in the descriptor form. AB13DX to compute the maximum singular value of a given continuous-time or discrete-time transfer-function matrix, either standard or in the descriptor form. AG07BD : to compute the inverse (Ai-lambda*Ei,Bi,Ci,Di) of a given descriptor system (A-lambda*E,B,C,D). DE01PD : to compute the convolution or deconvolution of two real signals using the Hartley transform. delctg : void logical function for LAPACK routine DGGES. DG01OD : to compute the (scrambled) discrete Hartley transform of a real signal. MA02HD : to check if A = DIAG*I, where I is an M-by-N matrix with ones on the diagonal and zeros elsewhere. MB02CU : to bring the first blocks of a generator to proper form. (Extended version of MB02CX, with a higher BLAS3 fraction and a pivoting scheme for rank-deficient generators.) MB02CV : to apply the transformations created by the SLICOT Library routine MB02CU on other columns / rows of the generator. MB02FD : to compute the incomplete Cholesky (ICC) factor of a symmetric positive definite block Toeplitz matrix T, defined by either its first block row, or its first block column. MB02GD : to compute the Cholesky factor of a banded symmetric positive definite block Toeplitz matrix. MB02HD : to compute the Cholesky factor of the matrix T'*T, with T a banded block Toeplitz matrix of full rank. MB02ID : to solve overdetermined or underdetermined real linear systems involving a full rank block Toeplitz matrix. MB02JD : to compute a full QR factorization of a block Toeplitz matrix of full rank. MB02JX : to compute a low rank QR factorization with column pivoting of a block Toeplitz matrix. MB02KD : to compute the product C = alpha*op( T )*B + beta*C, where alpha and beta are scalars and T is a block Toeplitz matrix. TG01BD : to reduce the matrices A and E of a system pencil corresponding to the descriptor triple (A-lambda E,B,C) to generalized upper Hessenberg form using orthogonal transformations. TG01WD : to reduce the pair (A,E) of a descriptor system to a real generalized Schur form using an orthogonal equivalence transformation, (A,E) <-- (Q'*A*Z,Q'*E*Z) and to apply the transformation to the matrices B and C: B <-- Q'*B and C <-- C*Z. New test programs: ================= TAB09JD : test program for AB09JD.f. TAB13DD : test program for AB13DD.f. TDE01PD : test program for DE01PD.f. TDG01OD : test program for DG01OD.f. TMB02FD : test program for MB02FD.f. TMB02GD : test program for MB02GD.f. TMB02HD : test program for MB02HD.f. TMB02ID : test program for MB02ID.f. TMB02JD : test program for MB02JD.f. TMB02JX : test program for MB02JX.f. TMB02KD : test program for MB02KD.f. Updated files for results: ========================= Using Compaq Visual Fortran 6.5, the following 4 .res files obtained on a PC machine, which previously slightly differ from those obtained on the Sun Solaris platform, now coincide (check!). These files are: AB09DD.res, BB03AD.res, SB08MD.res, and SB08ND.res. Documentation: ============= * : all newly added routines. libindex, : added links to the new document files. support New mexfiles: ============= bstred.f : balanced stochastic truncation model reduction based on SLICOT user-callable routine AB09HD. conred.f : frequency-weighted balancing related controller reduction based on SLICOT user-callable routine SB16AD. fwehna.f : frequency-weighted Hankel-norm approximation based on SLICOT user-callable routine AB09JD. fwered.f : frequency-weighted balancing related model reduction based on SLICOT user-callable routine AB09ID. sfored.f : state-feedback-observer controller reduction based on SLICOT user-callable routines SB16BD and SB16CD. linorm.f : L-infinity norm of a linear continuous- or discrete-time system based on SLICOT user-callable routine AB13DD. New m-files: ============ bst.m : balanced stochastic truncation based model reduction. cfconred.m : coprime factorization based controller reduction of state-feedback controllers. fwbconred.m : frequency-weighted balancing related controller reduction. fwbred.m : frequency-weighted balancing related model reduction. fwhna.m : frequency-weighted Hankel-notm approximation model reduction. sysredset.m : creation of a system order reduction options structure. sysredget.m : extracts the value of the named parameter from a system order reduction options structure. linorm.m : documentation for linorm.f mex-function. slinorm.m : L-infinity norm of a state-space system. test_conred.m : test files for various mexfiles and m-files. test_sysred_tools.m test_linorm.m test_slinorm.m conred.mat ======================================================================== Update February 25, 2001 March 3, 2001 Corrected Bugs: ============== AB09ED,: few lines added to finish if only unstable part is present. AB09MD, AB09ND AB13CD : adequate workspace has been provided even for the quick return cases; replaced IW3 by IW2 in the loop 30; replaced LDWORK by LDWORK-IWRK in the call of DGEES after the label 200; the routine now also deliver the peak frequency for which the estimated infinity norm is achieved; moreover, it works for systems without dynamic part (the case with n = 0 did not cause now the same quick return as for m = 0 or p = 0); the maximal singular value of D is properly taken into account; the lower bound on the norm is not set to 0 before the loop 320. SB02CX : changed the tolerance to 100*EPS. IB01ND : removed the first use of the non-initialized variable MAXWRK. IB01PD : replaced TOLL, at its first use, by TOLL1. recomputed RANKM before calling MB02QY. Updated Example Programs: ======================== TAB13CD : LIWORK is N; the unused variables LDAK, LDBK, LDCK, and LDDK have been deleted; the peak frequency is also displayed. New Routines: ============ AB09HD : to compute a reduced order model (Ar,Br,Cr,Dr) for an original state-space representation (A,B,C,D) by using the stochastic balancing approach in conjunction with the square-root or the balancing-free square-root Balance & Truncate (B&T) or Singular Perturbation Approximation (SPA) model reduction methods for the ALPHA-stable part of the system. AB09HX : to compute a reduced order model (Ar,Br,Cr,Dr) for an original stable state-space representation (A,B,C,D) by using the stochastic balancing approach in conjunction with the square-root or the balancing-free square-root Balance & Truncate (B&T) or Singular Perturbation Approximation (SPA) model reduction methods. The state dynamics matrix A of the original system is an upper quasi-triangular matrix in real Schur canonical form and D must be full row rank. AB09HY : to compute the Cholesky factors Su and Ru of the controllability Grammian P = Su*Su' and observability Grammian Q = Ru'*Ru. AB09ID : to compute a reduced order model (Ar,Br,Cr,Dr) for an original state-space representation (A,B,C,D) by using the frequency weighted square-root or balancing-free square-root Balance & Truncate (B&T) or Singular Perturbation Approximation (SPA) model reduction methods. The algorithm tries to minimize the norm of the frequency-weighted error ||V*(G-Gr)*W||, where G and Gr are the transfer-function matrices of the original and reduced order models, respectively, and V and W are frequency-weighting transfer-function matrices without poles on the boundary of the stability domain. If G is unstable, only the ALPHA-stable part of G is reduced. AB09IX : to compute a reduced order model (Ar,Br,Cr,Dr) for an original state-space representation (A,B,C,D) by using the square-root or balancing-free square-root Balance & Truncate (B&T) or Singular Perturbation Approximation (SPA) model reduction methods. The computation of truncation matrices is based on the Cholesky factor of a controllability Grammian and the Cholesky factor of an observability Grammian. AB09IY : to compute the Cholesky factors of the frequency-weighted controllability and observability Grammians corresponding to a frequency-weighted model reduction problem, for stable transfer-function matrices of the system and the weightings, with state matrices in real Schur form. AB09KD : to compute a reduced order model (Ar,Br,Cr,Dr) for an original state-space representation (A,B,C,D) by using the frequency weighted optimal Hankel-norm approximation method in conjunction with square-root balancing for the ALPHA-stable part of the system. AB09KX : to construct a state-space representation (A,BS,CS,DS) of the stable projection of V*G*W or conj(V)*G*conj(W) from the state-space representations of (A,B,C,D), (AV,BV,CV,DV), and (AW,BW,CW,DW) of the transfer-function matrices G, V and W, respectively, for the state matrix A in a real Schur form. It is assumed that G is stable. AB13MD : to compute an upper bound on the structured singular value for a given square complex matrix and a given block structure of the uncertainty. MB01WD : to compute the matrix formula _ R = alpha*( op( A )'*op( T )'*op( T ) + op( T )'*op( T )*op( A ) ) + beta*R, (1) or _ R = alpha*( op( A )'*op( T )'*op( T )*op( A ) - op( T )'*op( T )) + beta*R, (2) _ where alpha and beta are scalars, R, and R are symmetric matrices, T is a triangular matrix, A is a general or Hessenberg matrix, and op( M ) = M or op( M ) = M'. MB01XD : to compute the matrix product U' * U or L * L', where U and L are upper and lower triangular matrices, respectively, stored in the corresponding upper or lower triangular part of the array A (block algorithm). MB01XY : to compute the matrix product U' * U or L * L', where U and L are upper and lower triangular matrices, respectively, stored in the corresponding upper or lower triangular part of the array A (unblocked algorithm). MB01YD : to perform the symmetric rank k operations C := alpha*op( A )*op( A )' + beta*C, where alpha and beta are scalars, C is an n-by-n symmetric matrix, op( A ) is an n-by-k matrix, and op( A ) = A or op( A ) = A'. The matrix A has l nonzero codiagonals, either upper or lower. MB01ZD : to compute the matrix product H := alpha*op( T )*H, or H := alpha*H*op( T ), where alpha is a scalar, H is an m-by-n upper or lower Hessenberg-like matrix (with l nonzero subdiagonals or superdiagonals, respectively), T is a unit, or non-unit, upper or lower triangular matrix, and op( T ) = T or op( T ) = T'. MB02CD : to compute the Cholesky factor and the generator and/or the Cholesky factor of the inverse of a symmetric positive definite block Toeplitz matrix, defined by either its first block row, or its first block column. MB02CX : to bring the first blocks of a generator in proper form. MB02CY : to apply the transformations created by the SLICOT Library routine MB02CX on other columns or rows of the generator. MB02DD : to update the Cholesky factor, the generator, and/or the Cholesky factor of the inverse of a symmetric positive definite block Toeplitz matrix, given the information from a previous factorization and additional blocks of its first block row, or its first block column. MB02ED : to solve a system of linear equations, T*X = B or X*T = B, with a symmetric positive definite block Toeplitz matrix T, defined either by its first block row or its first block column. SB10ID : to compute the matrices of the positive feedback controller for the continuous-time shaped plant. SB10JD : to convert a descriptor state-space system into regular state-space form. SB10KD : to compute the matrices of the positive feedback controller for the discrete-time shaped plant. SB16AD : to compute a reduced order controller (Acr,Bcr,Ccr,Dcr) for an original state-space controller representation (Ac,Bc,Cc,Dc) by using the frequency weighted square-root or balancing-free square-root Balance & Truncate (B&T) or Singular Perturbation Approximation (SPA) model reduction methods. SB16AY : to compute the Cholesky factors of the frequency-weighted controllability and observability Grammians corresponding to a frequency-weighted model reduction problem, with the controller state matrix in a block-diagonal real Schur form. SB16BD : to compute, for a given open-loop model (A,B,C,D), and for a given state feedback gain F and full observer gain G, such that A+B*F and A+G*C are stable, a reduced order controller model (Ac,Bc,Cc,Dc) using a coprime factorization based controller reduction approach. SB16CD : to compute, for a given open-loop model (A,B,C,D), and for a given state feedback gain F and full observer gain G, such that A+B*F and A+G*C are stable, a reduced order controller model (Ac,Bc,Cc,Dc) using a coprime factorization based controller reduction approach. For reduction of the coprime factors, a stability enforcing frequency-weighted model reduction is performed using either the square-root or the balancing-free square-root versions of the Balance & Truncate (B&T) model reduction method. SB16CY : to compute, for a given open-loop model (A,B,C,0), and for a given state feedback gain F and full observer gain G, such that A+B*F and A+G*C are stable, the Cholesky factors Su and Ru of a controllability Grammian, P = Su*Su', and an observability Grammian, Q = Ru'*Ru, corresponding to a frequency-weighted model reduction of the left or right coprime factors of the state-feedback controller. Updated Routines (after February 26, 2001): ================ AB13MD : the input/output argument X was moved before the output arguments, and some improvements have been performed in the code. SB10ID : minor changes in the code and documentation. SB10JD minor changes in the documenting comments. SB10KD : improvements have been performed in the code, the workspace has been reduced. Updated files for results: ========================= For reference purposes, all example programs have been relinked with the prebuilt libraries for LAPACK compiled for Sun Solaris, SunOS enterprise 5.7 Generic_106541-03 sun4u sparc SUNW, Ultra-4 (file lapack_solaris.tgz), and LAPACK compiled for Pentium Pro, Windows NT, Visual Fortran (file lapack_win32.zip), both available on netlib, at the address http://www.netlib.org/lapack/archives. As a result, the following 25 .res files have been replaced: AB09AD.res, AB09BD.res, AB09CD.res, AB09ED.res, AB09FD.res, AB09GD.res, AB09MD.res, AB09ND.res, AB13FD.res, MB03QD.res, MB04DY.res, MB05MD.res, SB01BD.res, SB04MD.res, SB04PD.res, SB04QD.res, SB08CD.res, SB08DD.res, SB08ED.res, SB08FD.res, SB10FD.res, TB01KD.res, TB01LD.res, TB01TD.res, TB01WD.res. The differences to the previously stored values, obtained on a Sun Ultra 2 machine, using LAPACK Release 2.0 (with additional Release 3.0 codes for the generalized eigenproblems), are mainly due to the different Schur forms computed. Namely, the 2-by-2 blocks could have the off-diagonal elements permuted, or with the signs interchanged, and/or the diagonal blocks could come in a different order, in which case the corresponding rows and columns have different elements. The differences are essentially those described at the previous update, in the file readme from the compressed tar archive diff_results.tar.gz, then stored in the root directory of the SLICOT ftp site: - AB09AD, AB09BD, AB09FD, AB09GD, SB04MD, TB01KD, TB01LD, and TB01WD : differ only in the signs of some elements; - AB13FD : one or two less significant digits differ; - AB09CD, AB09ED, AB09MD, AB09ND, and MB03QD : some rows/columns are permuted, possibly with signs of some elements changed; - SB01BD : a different relative residual norm (the absolute residual norm is close to machine accuracy); - MB04DY : different scaling factors (and, hence, off-diagonal elements); - TB01TD : different balancing factors; - SB08CD, SB08DD, SB08ED, and SB08FD : Schur forms are different, as mentioned above; - MB05MD : some eigenvalues are interchanged, and the corresponding rows/columns are permuted; - SB10FD : two real poles are interchanged; - SB04PD and SB04QD : two columns of the orthogonal matrix differ. The following 5 .res files obtained on a PC machine slightly differ from those obtained on the Sun Solaris platform: AB09DD.res, BB03AD.res, MC03ND.res, SB08MD.res, and SB08ND.res: - SB08MD and SB08ND : differ only in white spaces; - MC03ND : differs only in the signs of some elements; - BB03AD and AB09DD : one or two less significant digits differ. The file slicotPC.zip contains the results got on a PC machine. Documentation: ============= AB09ND : On exit, if INFO = 0, IWORK(1) contains the order of the minimal realization of the ALPHA-stable part of the system. AB13CD : inserted a function value description; IWORK has length just N; on exit, DWORK(2) now contains the frequency where the gain of the frequency response achieves its peak value. MB02SD,: changed the operations count to 0(n^2) at NUMERICAL ASPECTS. MB02SZ libindex, : added links to the new document files. support The documentation files corresponding to the changed .res files mentioned above have been also updated. The results on the Sun platform were used. Corrected Bugs in mexfiles: =========================== aresolc.f : an error is produced if the reciprocal condition number and forward error bound are asked, but cannot be computed; these results are obtained just for meth = 1, and flag(4) = 2; replaced the name ARESOL by ARESOLC. findBD.f : moved the tests involving IJOB after IJOB initialization. order.f : replaced > by .GT. in a test for CSIZE. sysred.f : replaced MAX(N, M, P) by MAX(N, M+P) in the formula for LWR. Corrected Bugs in m-files: ========================== test_aresol.m : replaced find(ev ... by ~isempty(find(real(ev) ... test_aresolc.m and find(abs(ev) ... by ~isempty(find(abs(ev) ... New mexfiles: ============= clsdp.f : Loop Shaping Design of continuous-time systems dlsdp.f : Loop Shaping Design of discrete-time systems fstoep.f : fast Toeplitz solver for symmetric positive definite matrices mucomp.f : computing the structured singular value New m-files: ============ clsdp.m, : m-files for the loop shaping design dlsdp.m, (added after February 26, 2001) tclsdp.m, tdlsdp.m fstoep.m, : m-files for the fast Toeplitz solver. fstchol.m, fstgen.m, fstsol.m, fstupd.m, test_fstoep.m, test_fstoem.m mucomp.m, : m-files for the structured singular value problem tmucomp.m (added after February 26, 2001) Updated Mexfiles: ================= hinorm : added the peak frequency for which the estimated infinity norm is achieved as an optional output parameter; TOL cannot be smaller than 100*EPS, where EPS is the relative machine precision; removed AB13CD from the External Subroutines list and added AB13CD and DLAMCH in the External Functions list. New Unix mex- and m-files: (after February 26, 2001) ========================== The previous Unix mexfiles used maximally dimensioned arrays, according to the Fortran 77 standards. The previous archives in the subdirectories MatlabTools/Unix/SLmex and MatlabTools/Unix/SLTools have been renamed as mexfiles77.tar.gz and mfiles77.tar.gz, respectively. They will not be updated in the future. New archives, mexfiles.tar.gz and mfiles.tar.gz, which use the Fortran 90/95 dynamic allocation feature for arrays, have been added in the same subdirectories. These archives have now essentially the same contents as the corresponding Windows archives, mexfiles.zip and mfiles.zip, respectively. The compressed Unix library file slicot.tar.gz includes now subdirectories called F77, added to both SLmex and SLTools subdirectories. The Fortran 77-related files have been moved into these F77 subdirectories, while the SLmex and SLTools subdirectories contain the Fortran 90/95 versions of the mexfiles and associated m-files, respectively. Moreover, the mexfiles have been also compiled for Matlab 6 (using Compaq Visual Fortran Version 6.5), and they are stored in the archive dllfiles6.zip, subirectory MatlabTools/Windows/SLTools. New contributed m-files (included in the LYAPACK collections in the ======================= subdirectory MatlabTools/contrib) lyapack.tar.gz : LYAPACK collection for Unix platforms (see below); lyapack.zip : LYAPACK collection for Windows platforms; lyapack.pdf : Adobe Acrobat pdf Users' Guide documentation file for the LYAPACK collection. Note: The two compressed files (.gz and .zip) contain the Users' Guide in both Postscript and pdf formats. LYAPACK is an efficient MATLAB toolbox (based on MATLAB m-functions) for the solution of certain large scale problems in control theory, which are closely related to Lyapunov equations. It can solve Lyapunov and Riccati equations, and do model reduction. LYAPACK uses iterative algorithms and is intended for large, sparse problems. LYAPACK Developer: Dr. Thilo Penzl. New demo files: ============== fstoep_demo.zip, : Matlab 5.3 demonstration packages for the SLICOT slident_demo.zip fast Toeplitz solvers and linear system identification software. ======================================================================== Update September 21, 2000 Corrected Bugs: ============== AB07ND : replaced ONE by -ONE in the call of DGEMV in the loop 20. IB01AD,: added (M+L)*2*NOBR in the expression for LDWORK when ALG = 'F' IB01MD and BATCH = 'O' or BATCH = 'L' and CONCT = 'N'; reduced the length of LDWORK in the other cases for ALG = 'F'. IB01MD : activated the call of IB01MY (for ALG = 'F'); return after the call of IB01MY if BATCH = 'F' or 'I'. SB04MD : set DWORK(1) = ONE for Quick return, and checked that LDWORK is at least 1. SB04QD : checked that LDWORK is at least 1. Updated Routines: ================ IB01OD : replaced the second parameter in the call of IB01OY with NOBR-1. Also, the estimated order is not allowed to exceed NOBR-1. Updated Example Programs: ======================== TAB01ND,: added the type declaration "INTRINSIC MAX". TAB01OD, TAB05RD, TAB07MD, TAB08ND, TAB09ED, TAB09FD, TAB09GD, TAB09MD, TAB09ND, TAB13AD, TAB13BD, TAB13CD, TFB01QD, TFB01RD, TFB01SD, TFB01TD, TFB01VD, TIB01AD, TIB01BD, TIB01CD, TMB02MD, TMB02QD, TMB03UD, TMB04OD, TMB04UD, TMB04VD, TMB04XD, TMC01RD, TMC03ND, TSB01BD, TSB01DD, TSB02ND, TSB02OD, TSB02PD, TSB02RD, TSB03OD, TSB03QD, TSB03SD, TSB03TD, TSB03UD, TSB04MD, TSB04ND, TSB04OD, TSB04PD, TSB04QD, TSB04RD, TSB06ND, TSB08CD, TSB08DD, TSB08ED, TSB08FD, TSB10DD, TSB10ED, TSB10FD, TSB10HD, TSG03AD, TSG03BD, TTB01PD, TTB01UD, TTB01ZD, TTC01OD, TTG01AD, TTG01CD, TTG01DD, TTG01ED, TTG01FD, TTG01HD, TTG01ID, TTG01JD TAB09ED,: added the type declaration "INTRINSIC MIN". TAB13BD, TIB01CD, TMB02QD, TMB03OD, TMB03PD, TMB04GD, TSB03OD, TTG01CD, TTG01DD, TTG01ED, TTG01FD, Note: the changes related to the INTRINSIC declarations above are not incorporated in the module files AB01ND.tar.gz, etc., except when additional changes were needed. TSB04MD,: LDWORK is at least 1. TSB04QD New Routines: ============ IB01MY : to compute the upper triangular factor in the QR factorization of the concatenated block Hankel matrices, using a fast QR algorithm. MA02FD : to construct a hyperbolic plane rotation. New Directory (under the subdirectory MatlabTools) ============= contrib : contains externally contributed Matlab files and mexfiles based on SLICOT library. Included are, currently, a Fortran program, two routines, and other associated files for calling SLICOT TB03AD routine from a Matlab environment. These files were developed by Dr. P.-O. Malaterre and used to compare SLICOT results with the Polynomial Toolbox results. The codes are described in the NICONET Newsletter, 5, pp. 13-20, July 2000. Documentation: ============= IB01AD,: corrected the length of DWORK and its description for ALG = 'F'; IB01MD corrected the factor K (for ALG = 'F') of the number of connection elements (used in the description of DWORK). IB01BD : the length of IWORK is now also specified for METH = 'C'. SB04MD,: LDWORK is at least 1. SB04QD libindex, : added links to the new document files. support Corrected Bugs in m-files: ========================== markov.m : the previously omitted term C*B was added, and the revised code is more efficient. findABCD.m, : few comments have been changed. findAC.m, findBD.m, findBDK.m, findR.m, inistate.m, order.m, sident.m, slgely.m, slgesg.m, slgest.m, slgsly.m, slgsst.m, slsylv.m findx0BD.m : few comments have been changed, and the tests for withx0 and withd have been corrected. inistate.m : the value for job when calling findBD for nargin = 6 has been corrected. Improved m-files: ================ slsylv.m : a new option enables to choose between Bartels-Stewart and Hessenberg-Schur methods. New m-files: ============ sldsyl.m : to solve discrete-time Sylvester equations. slmoesp.m : to identify a linear system using MOESP approach. sln4sid.m : to identify a linear system using N4SID approach. slmoen4.m : to identify a linear system using combined MOESP and N4SID approaches. slmoesm.m : to identify a linear system using combined MOESP and simulation approaches. test_order : tests order.f for multiple data batches. call_torder: calls test_order with different option combinations. test_sl* : similarly, for * = moen4, moesm, moesp, or n4sid. call_tsl* Updated Mexfiles: ================= linmeq : discrete-time Sylvester equations can now be solved; Sylvester equations with op(A) = 'N', op(B) = 'T', or op(A) = 'T', op(B) = 'N', can be solved. A new option for Sylvester equations enables to select between the Bartels-Stewart and Hessenberg-Schur methods. order : few changes have been performed mainly for computing LDWORK and LIWORK, and for saving the connection elements. sident : few changes in the comments, and for computing LDWORK for METH = 'C'. findBD : few changes in the comments, as well as for returning X0. Also, the first dimension of U is not checked when M = 0. Alternate Results: ================= Recently, it has been drawn to our attention that some results obtained using the SLICOT example programs could (slightly) depend on the specific platform and LAPACK release installed. The files included in the standard SLICOT collection (the compressed archives slicot.tar.gz, slicotPC.zip, AB01MD.tar.gz, etc., as well as in the documentation html files) have been obtained on a Sun Ultra 2 machine, using LAPACK Release 2.0 (with additional Release 3.0 codes for the generalized eigenproblems). The differences that have been observed when using a PC and LAPACK Release 3.0 are mainly due to the different Schur forms computed. Namely, the 2-by-2 blocks could have the off-diagonal elements permuted, or with the signs interchanged, and/or the diagonal blocks could come in a different order, in which case the corresponding rows and columns have different elements. For convenience, the alternate results obtained on a PC with LAPACK Release 3.0 are stored in the compressed tar file diff_results.tar.gz, in the root directory of the SLICOT ftp site. Details are given in the readme file included in the mentioned archive. ======================================================================== Update June 30, 2000 Corrected Bugs: ============== AB09BX : if the order of the reduced system is zero, call AB09DD and exit. AB09CX : the system is transformed such that Ar is in a real Schur form also for NR.EQ.NMINR. MB01TD : the non-initialized value of DWORK(N) is no longer used. The length of DWORK is N-1. SB04MD : replaced 3*M by 5*M in the workspace length. SB04ND SB08HD : replaced the call of LAPACK routine DLAPMT by a call to the new SB10HD SLICOT routine MA02GD. TB03AD : replaced the call of LAPACK routine DLAPMT by a call to the new SLICOT routine MA02GD; corrected the description of the parameter QCOEFF when LERI ='R'; TTB03AD has been updated accordingly. Improved Routines: ================= SB03MX : replaced the calls to SLICOT routine SB03MU by either calls to LAPACK routine DLALN2 (for linear systems of order 2), or to the new, slightly more efficient SLICOT routine SB04PX (for linear systems of order 4). SB03OR : replaced the called routine SB03MU by the new, slightly more efficient routine SB04PX. SB04ND : calls DTRSYL when both matrices A and B are in real Schur form. Updated Routines, due to the changes in MB03UD and in the LAPACK routines ================ DBDSQR and DGESVD (in LAPACK Release 3.0, Update October 31, 1999): AB09AX : replaced MAX(4*N, 5*N-4) by 5*N in a comment concerning AB09BX, MB03UD. AB13AX MB02MD : replaced everywhere (in the workspace statements) "5*x-4" by "5*x", due to the change in DGESVD. MB02UD : replaced everywhere "4*L, 5*L-4" by "5*L", due to the change in MB03UD. MB03UD : replaced LDWORK = MAX(1, 4*N, 5*N-4) by LDWORK = MAX(1, 5*N), due to the change in DBDSQR. SB10DD, : replaced everywhere (in the workspace statements) "5*x-4" by SB10ED, "5*x", due to the change in DGESVD. SB10FD, SB10HD, SB10PD, SB10UD TG01ED : replaced 5*MIN(L,N)-4 by 5*MIN(L,N) in LDWORK, due to the changes in MB03UD and in DGESVD. Updated Example Programs: ======================== TMB01TD, TMB02MD : updated the length of LDWORK, due to the changes in TMB03UD, TSB04MD MB01TD, MB02MD, MB03UD, SB04MD, SB04ND, and TG01ED, TSB04ND, TTG01ED respectively. TTB03AD : corrected printing of QCOEFF when LERI ='R'. New Routines: ============ AB07ND : compute the inverse (Ai,Bi,Ci,Di) of a given system (A,B,C,D). IB01AD : preprocess the input-output data for estimating the matrices of a linear time-invariant dynamical system and find an estimate of the system order. The input-output data can, optionally, be processed sequentially. IB01BD : estimate the system matrices A, C, B, and D, the noise covariance matrices Q, Ry, and S, and the Kalman gain matrix K of a linear time-invariant state space model, using the processed triangular factor R of the concatenated block-Hankel matrices, provided by SLICOT Library routine IB01AD. IB01CD : estimate the initial state and, optionally, the system matrices B and D of a linear time-invariant discrete-time system, given the system matrices (A,B,C,D), or only the matrix pair (A,C), and the input and output trajectories of the system. IB01MD : compute the upper triangular factor in the QR factorization of the block-Hankel matrix built from input-output data for subspace identification. IB01ND : find the singular value decomposition giving the system order, using the triangular factor of the concatenated block-Hankel matrices. IB01OD : find the system order, using the singular value decomposition. IB01OY : ask for user's confirmation of the system order found. IB01PD : estimate the system matrices and covariance matrices of a linear time-invariant state space model. IB01PX : compute the matrices B and D of a linear time-invariant state space model using Kronecker products. IB01PY : compute the triangular QR factor of a structured matrix and estimate the matrices B and D of a linear time-invariant state space model. IB01QD : estimate the initial state and the system matrices B and D of a linear time-invariant state space model, given the matrix pair (A,C) and the input and output trajectories. IB01RD : estimate the initial state of a linear time-invariant state space model, given the system matrices (A,B,C,D) and the input and output trajectories. MA02GD : perform a series of column interchanges on the matrix factored by LAPACK Library routine DGETRF. MB01VD : compute the Kronecker product of two matrices. MB02VD : solve X*A = B or X*A' = B. SB02PD : solve continuous-time Riccati equations by the matrix sign function and estimate the reciprocal condition number and forward error bound. SB04PD : solve continuous-time or discrete-time Sylvester equations. SB04PX : solve discrete-time Sylvester equations with matrices of order at most 2 (slightly improved version of SB03MU). SB04PY : solve discrete-time Sylvester equations with matrices in real Schur form. SB04QD : solve discrete-time Sylvester equations using the Hessenberg-Schur method. SB04QR : solve a linear algebraic system of order M whose coefficient matrix has zeros below the third subdiagonal and zero elements on the third subdiagonal with even column indices. The matrix is stored compactly, row-wise. SB04QU : construct and solve a linear algebraic system of order 2*M whose coefficient matrix has zeros below the third subdiagonal, and zero elements on the third subdiagonal with even column indices. Such systems appear when solving discrete-time Sylvester equations using the Hessenberg-Schur method. SB04QY : construct and solve a linear algebraic system of order M whose coefficient matrix is in upper Hessenberg form. Such systems appear when solving discrete-time Sylvester equations using the Hessenberg-Schur method. SB04RD : solve for X the discrete-time Sylvester equation X + AXB = C, with at least one of the matrices A or B in Schur form and the other in Hessenberg or Schur form (both either upper or lower). SB04RV : construct the right-hand sides for a system of equations in Hessenberg form solved via SB04RX (case with 2 right-hand sides). SB04RW : construct the right-hand side for a system of equations in Hessenberg form solved via SB04RY (case with 1 right-hand side). SB04RX : solve a system of equations in Hessenberg form with two consecutive offdiagonals and two right-hand sides (discrete-time case). SB04RY : solve a system of equations in Hessenberg form with one offdiagonal and one right-hand side (discrete-time case). Documentation: ============= MB01TD, MB02MD, : updated the length of LDWORK, due to the changes in MB02UD, MB03UD, corresponding subroutines. SB10DD, SB10ED, SB10FD, SB10HD, SB10PD, SB10UD, TG01ED SB04ND : TOL, IWORK and DWORK are not needed when both matrices A and B are in real Schur form. TB03AD : corrected the description of the parameter QCOEFF when LERI ='R'. New Mexfiles: ============ Three new mexfiles and ten Matlab m files have been added in the subdirectories ./SLmex and ./SLTools, respectively. They are useful for performing system identification using subspace techniques (MOESP, N4SID, or their combination). ======================================================================== Update February 18, 2000 Known Bugs: =========== TB03AD: Few days ago, a bug has been reported. It has been almost fixed, but additional analysis and improvements are needed. The routine will be replaced by the corrected one as soon as possible. Bugs Corrected: ============== TAB13ED: deleted the second definition of the variable LDA. TAB13FD TAG08BD: changed a format statement to avoid a warning when using F90 TMB03WD compiler on a SUN machine. New Routines Added: ================== BB01AD : generate the benchmark examples for the numerical solution of continuous-time algebraic Riccati equations (CAREs). The associated data files have been added in the subdirectory benchmark_data. BB02AD : generate the benchmark examples for the numerical solution of discrete-time algebraic Riccati equations (DAREs). The associated data files have been added in the subdirectory benchmark_data. New Data Files Added: ==================== BB01*.dat: 14 data files to be used in conjunction with the new Riccati BB02*.dat benchmark collections have been added in the subdirectory benchmark_data. Documentation: ============= SB10HD : deleted the compound matrices in the description of the error return INFO = 3. SB10TD : replaced D2 by D22 in the formula from the NUMERICAL ASPECTS section. libindex : updated in the files slicot.tar.gz and slicotPC.zip. New Mexfiles Added: =================== aresolc.f: added (in conjunction with aresolc.m, test_aresolc.m). The help file for the standard Riccati solver (aresol.m) has been modified. New Makefiles Added: =================== Makefiles for compiling the library, as well as for linking and running the example programs on PC's under Windows 95, 98, or NT operating systems, have been added. ======================================================================== Update December 18, 1999 Bugs Corrected: ============== SB04OD : initialized IJOB = 0 also when TRANS = 'T', to avoid an error message for a non-initialized variable on some compilers. TAB09CD: deleted the last " )" in "( NMAX*( NMAX + 1 ) )/2 )" from the statement defining LDWORK. TMB03OD: replaced "Row" by "Columns" in FORMAT statement 99994. New Routines Added: ================== AB13ED : estimates the distance from a real matrix to the nearest complex matrix with an eigenvalue on the imaginary axis, using bisection. AB13FD : computes the distance from a real matrix to the nearest complex matrix with an eigenvalue on the imaginary axis, using bisection and SVD (quadratically convergent algorithm). FD01AD : computes a fast recursive least-squares filter using a QR-decomposition based approach. MB03NY : computes the smallest singular value of A - jwI. Changes in the organization of the compressed library files: =========================================================== The directory tree structure of the compressed library files slicot.tar.gz, for Unix platforms, and slicotPC.zip (former slicotPC.tar.gz), for Windows platforms (Windows 95, 98, and NT), has been changed. The new tree structure is as follows: slicot : contains the files libindex.html, make.inc, makefile, and the following subdirectories: benchmark_data : contains benchmark data files for Fortran benchmark routines (*.dat); doc : contains SLICOT documentation files for routines, (.html); examples : contains SLICOT example programs, data, and results (.f, .dat, .res), and makefile, for compiling, linking and executing these programs; src : contains SLICOT source files for routines (.f), and makefile, for compiling all routines; SLTools : contains Matlab .m files, mex files (optionally), and data .mat files (optionally); SLmex : contains Fortran source codes for mexfiles. Currently, there are no makefiles for the PC version, and no executable Matlab files (e.g., .dll files) have been included. These files may, however, be separately taken from the ftp site. New Directories Added: ===================== FD/ : Fast Recursive Least Squares Filters. FD01/ MatlabTools/ : Matlab files and mexfiles, contained in its subdirectories: Unix/ : Matlab files and mexfiles for Unix platforms. SLTools/ : Matlab source and data files (.m and .mat). SLmex/ : Fortran source codes for Matlab mexfiles based on SLICOT library. Windows/ : Matlab files and mexfiles for Windows platforms. SLToolboxes/ : Matlab executable files (.dll) groupped in several archives, according to their functionality. SLTools/ : Matlab source and data files (.m and .mat). SLdemos/ : Matlab source and executable files (.m, .mat, and .dll), for a GUI-based demo. SLmex/ : Fortran source codes for Matlab mexfiles based on SLICOT library. contrib/ : routines which were previously proposed for being included in the SLICOT Library, but which could not be included yet. These routines do not follow the latest SLICOT Implementation and Documentation Standards, but could be of interest to some users. Included are G. Miminis' routines for solving the Pole Placement Single-input or Multi-input problem, and A.J. Geurts and C. Praagman's routines for computing a unimodular polynomial matrix U(s) such that R(s) = P(s) * U(s) is column reduced, given a polynomial matrix P(s). Note: the former directory Mexfiles (under SLICOT/ directory), and all its subdirectories have been removed, their contents being replaced by the contents of the new directory MatlabTools/ and its subdirectories. ======================================================================== Update November 27, 1999 Important Notice: Release 4.0 of SLICOT Library took place =========================================================== Release 4.0 of the SLICOT Library has been announced and took place in November 27, 1999. All Library source files then contained the Release 4.0 statement notice. Several routines, together to the associated documentation files, example program files, data and results, have been removed from SLICOT Release 3.0. The removed files have better counterparts, and their removal has been announced in the previous versions of the Release.Notes file (see Release3.History file). The table below indicates the removed files, and their functionally equivalent routines. ----------------------------------------------------------------------- Table 1: Routines removed in Release 4.0 of the SLICOT Library and their functional correspondents. ----------------------------------------------------------------------- Routines removed Functionally equivalent routines ----------------------------------------------------------------------- AB06MD.f TB01PD.f MB04SD.f MB04UD.f MB04TD.f MB04VD.f TB01QD.f TB04AD.f TB01QY.f TB04AY.f TB01RD.f TB05AD.f TB01SD.f TB03AD.f TB01SY.f TB03AY.f TC01MD.f TC05AD.f TC01ND.f TC04AD.f TD01MD.f TD05AD.f TD01ND.f TD03AD.f TD01NY.f TD03AY.f TD01OD.f TD04AD.f ----------------------------------------------------------------------- In addition, the entire subdirectory LAPACK_N of /pub/WGS/SLICOT/, containing used LAPACK files not included in LAPACK Release 2.0, has been removed, because the needed files are now included in the LAPACK Release 3.0. ========================================================================