Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[8752] trunk/octave-forge/extra/control-devel
|
From: <par...@us...> - 2011-10-16 09:47:10
|
Revision: 8752
http://octave.svn.sourceforge.net/octave/?rev=8752&view=rev
Author: paramaniac
Date: 2011-10-16 09:47:01 +0000 (Sun, 16 Oct 2011)
Log Message:
-----------
control: add identification routines
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/src/MA02AD.f
trunk/octave-forge/extra/control-devel/src/MB01SD.f
trunk/octave-forge/extra/control-devel/src/MB02UD.f
trunk/octave-forge/extra/control-devel/src/MB03UD.f
trunk/octave-forge/extra/control-devel/src/MB04OD.f
trunk/octave-forge/extra/control-devel/src/MB04OY.f
trunk/octave-forge/extra/control-devel/src/SB02MR.f
trunk/octave-forge/extra/control-devel/src/SB02MS.f
trunk/octave-forge/extra/control-devel/src/SB02MV.f
trunk/octave-forge/extra/control-devel/src/SB02MW.f
trunk/octave-forge/extra/control-devel/src/SB04PX.f
trunk/octave-forge/extra/control-devel/src/TB01WD.f
trunk/octave-forge/extra/control-devel/src/select.f
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/DESCRIPTION
trunk/octave-forge/extra/control-devel/INDEX
trunk/octave-forge/extra/control-devel/devel/makefile_conred.m
trunk/octave-forge/extra/control-devel/devel/makefile_ident.m
trunk/octave-forge/extra/control-devel/devel/makefile_modred.m
trunk/octave-forge/extra/control-devel/src/IB01AD.f
trunk/octave-forge/extra/control-devel/src/IB01BD.f
trunk/octave-forge/extra/control-devel/src/IB01CD.f
trunk/octave-forge/extra/control-devel/src/IB01MD.f
trunk/octave-forge/extra/control-devel/src/IB01MY.f
trunk/octave-forge/extra/control-devel/src/IB01ND.f
trunk/octave-forge/extra/control-devel/src/IB01OD.f
trunk/octave-forge/extra/control-devel/src/IB01OY.f
trunk/octave-forge/extra/control-devel/src/IB01PD.f
trunk/octave-forge/extra/control-devel/src/IB01PX.f
trunk/octave-forge/extra/control-devel/src/IB01PY.f
trunk/octave-forge/extra/control-devel/src/IB01QD.f
trunk/octave-forge/extra/control-devel/src/IB01RD.f
trunk/octave-forge/extra/control-devel/src/MA02ED.f
trunk/octave-forge/extra/control-devel/src/MA02FD.f
trunk/octave-forge/extra/control-devel/src/MB01RU.f
trunk/octave-forge/extra/control-devel/src/MB01RX.f
trunk/octave-forge/extra/control-devel/src/MB01RY.f
trunk/octave-forge/extra/control-devel/src/MB01TD.f
trunk/octave-forge/extra/control-devel/src/MB01UD.f
trunk/octave-forge/extra/control-devel/src/MB01VD.f
trunk/octave-forge/extra/control-devel/src/MB02PD.f
trunk/octave-forge/extra/control-devel/src/MB02QY.f
trunk/octave-forge/extra/control-devel/src/MB03OD.f
trunk/octave-forge/extra/control-devel/src/MB04ID.f
trunk/octave-forge/extra/control-devel/src/MB04IY.f
trunk/octave-forge/extra/control-devel/src/MB04KD.f
trunk/octave-forge/extra/control-devel/src/Makefile
trunk/octave-forge/extra/control-devel/src/SB02MT.f
trunk/octave-forge/extra/control-devel/src/SB02ND.f
trunk/octave-forge/extra/control-devel/src/SB02QD.f
trunk/octave-forge/extra/control-devel/src/SB02RD.f
trunk/octave-forge/extra/control-devel/src/SB02RU.f
trunk/octave-forge/extra/control-devel/src/SB02SD.f
trunk/octave-forge/extra/control-devel/src/SB03MV.f
trunk/octave-forge/extra/control-devel/src/SB03MW.f
trunk/octave-forge/extra/control-devel/src/SB03MX.f
trunk/octave-forge/extra/control-devel/src/SB03MY.f
trunk/octave-forge/extra/control-devel/src/SB03QX.f
trunk/octave-forge/extra/control-devel/src/SB03QY.f
trunk/octave-forge/extra/control-devel/src/SB03SX.f
trunk/octave-forge/extra/control-devel/src/SB03SY.f
Removed Paths:
-------------
trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m
trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m
Added: trunk/octave-forge/extra/control-devel/DESCRIPTION
===================================================================
--- trunk/octave-forge/extra/control-devel/DESCRIPTION (rev 0)
+++ trunk/octave-forge/extra/control-devel/DESCRIPTION 2011-10-16 09:47:01 UTC (rev 8752)
@@ -0,0 +1,11 @@
+Name: Control-Devel
+Version: 0.1.50
+Date: 2011-10-16
+Author: Lukas Reichlin <luk...@gm...>
+Maintainer: Lukas Reichlin <luk...@gm...>
+Title: Control Systems Developer's Playground
+Description: SLICOT system identification plus model and controller reduction
+Depends: octave (>= 3.3.90), control (>= 2.2.0)
+Autoload: yes
+License: GPL version 3 or later
+Url: http://octave.sf.net, http://www.slicot.org
Added: trunk/octave-forge/extra/control-devel/INDEX
===================================================================
--- trunk/octave-forge/extra/control-devel/INDEX (rev 0)
+++ trunk/octave-forge/extra/control-devel/INDEX 2011-10-16 09:47:01 UTC (rev 8752)
@@ -0,0 +1,14 @@
+control-devel >> Control Theory
+Examples
+Experimental Data Handling
+ @iddata/get
+ iddata
+ @iddata/set
+ @iddata/size
+System Identification
+ moesp
+ n4sid
+Model Reduction
+ hna
+Controller Reduction
+
Copied: trunk/octave-forge/extra/control-devel/devel/makefile_conred.m (from rev 8751, trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m)
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/makefile_conred.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/makefile_conred.m 2011-10-16 09:47:01 UTC (rev 8752)
@@ -0,0 +1,26 @@
+homedir = pwd ();
+develdir = fileparts (which ("makefile_conred"));
+srcdir = [develdir, "/../src"];
+cd (srcdir);
+
+mkoctfile SB16AD.f TB01ID.f SB16AY.f TB01KD.f AB09IX.f \
+ MB04OD.f MB01WD.f SB03OD.f MB03UD.f AB05PD.f \
+ AB09DD.f AB07ND.f TB01LD.f AB05QD.f SB03OU.f \
+ MA02AD.f MB03QX.f select.f MB01YD.f MB01ZD.f \
+ SB03OT.f MB04OY.f MB03QD.f MB04ND.f MB03QY.f \
+ SB03OR.f SB03OY.f SB04PX.f MB04NY.f SB03OV.f
+
+mkoctfile SB16BD.f AB09AD.f AB09BD.f SB08GD.f SB08HD.f \
+ TB01ID.f AB09AX.f MA02GD.f AB09BX.f TB01WD.f \
+ MA02DD.f MB03UD.f select.f AB09DD.f SB03OU.f \
+ MA02AD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \
+ SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f
+
+mkoctfile SB16CD.f SB16CY.f AB09IX.f SB03OD.f MB02UD.f \
+ AB09DD.f MA02AD.f MB03UD.f select.f SB03OU.f \
+ MB01SD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \
+ SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f
+
+system ("rm *.o");
+cd (homedir);
+
Added: trunk/octave-forge/extra/control-devel/devel/makefile_ident.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2011-10-16 09:47:01 UTC (rev 8752)
@@ -0,0 +1,28 @@
+homedir = pwd ();
+develdir = fileparts (which ("makefile_ident"));
+srcdir = [develdir, "/../src"];
+cd (srcdir);
+
+## preprocess the input-output data
+mkoctfile IB01AD.f IB01MD.f IB01ND.f IB01OD.f IB01MY.f \
+ MB04OD.f MB03UD.f MB04ID.f MA02AD.f MB03OD.f \
+ MB04IY.f IB01OY.f MA02ED.f MA02FD.f MB04OY.f
+
+## estimating system matrices, Kalman gain, and covariances
+mkoctfile IB01BD.f IB01PD.f MA02AD.f SB02MT.f SB02RD.f \
+ SB02ND.f MB02UD.f MA02ED.f IB01PY.f MB03OD.f \
+ MB02QY.f IB01PX.f SB02MS.f SB02RU.f SB02SD.f \
+ MB01RU.f SB02QD.f SB02MV.f SB02MW.f SB02MR.f \
+ MB02PD.f MB01SD.f MB04KD.f MB03UD.f MB04OD.f \
+ MB04OY.f MB01VD.f select.f MB01UD.f SB03SY.f \
+ MB01RX.f SB03MX.f SB03SX.f MB01RY.f SB03QY.f \
+ SB03QX.f SB03MY.f SB04PX.f SB03MV.f SB03MW.f
+
+## estimating the initial state
+mkoctfile IB01CD.f TB01WD.f IB01RD.f IB01QD.f select.f \
+ MB01TD.f MA02AD.f MB04OD.f MB04OY.f MB02UD.f \
+ MB03UD.f MB01SD.f
+
+system ("rm *.o");
+cd (homedir);
+
Copied: trunk/octave-forge/extra/control-devel/devel/makefile_modred.m (from rev 8751, trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m)
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/makefile_modred.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/makefile_modred.m 2011-10-16 09:47:01 UTC (rev 8752)
@@ -0,0 +1,35 @@
+homedir = pwd ();
+develdir = fileparts (which ("makefile_modred"));
+srcdir = [develdir, "/../src"];
+cd (srcdir);
+
+mkoctfile AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \
+ AB09IX.f AB09IY.f SB08CD.f MB04ND.f TB01XD.f \
+ MB04OD.f MB01WD.f MB03UD.f AB07MD.f SB01FY.f \
+ AB09DD.f TB01LD.f SB03OU.f TB01UD.f MA02AD.f \
+ MA02BD.f MB03OY.f MB03QX.f MB01PD.f select.f \
+ MB01YD.f MB04NY.f MB01ZD.f SB03OT.f MB04OX.f \
+ MB04OY.f MB03QD.f SB03OY.f MB03QY.f MB01QD.f \
+ SB03OR.f SB03OV.f SB04PX.f
+
+mkoctfile "-Wl,-framework" "-Wl,vecLib" \
+ slab09jd.cc \
+ AB09JD.f TB01ID.f TB01KD.f AB07ND.f AB09JV.f \
+ AB09JW.f AB09CX.f AG07BD.f AB08MD.f AB04MD.f \
+ TB01LD.f delctg.f SB04PY.f AB09AX.f AB08NX.f \
+ MB01SD.f AB09JX.f MA02AD.f TB01WD.f MB03OY.f \
+ MB03PY.f MA02DD.f MB03UD.f MB03QX.f select.f \
+ SB04PX.f SB03OU.f MB03QD.f MB03QY.f SB03OT.f \
+ MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \
+ MB04OY.f SB03OV.f
+
+mkoctfile AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \
+ AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \
+ SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \
+ SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \
+ SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \
+ SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f
+
+system ("rm *.o");
+cd (homedir);
+
Deleted: trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m 2011-10-15 19:22:25 UTC (rev 8751)
+++ trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m 2011-10-16 09:47:01 UTC (rev 8752)
@@ -1,26 +0,0 @@
-homedir = pwd ();
-develdir = fileparts (which ("makefile_slconred"));
-srcdir = [develdir, "/../src"];
-cd (srcdir);
-
-mkoctfile SB16AD.f TB01ID.f SB16AY.f TB01KD.f AB09IX.f \
- MB04OD.f MB01WD.f SB03OD.f MB03UD.f AB05PD.f \
- AB09DD.f AB07ND.f TB01LD.f AB05QD.f SB03OU.f \
- MA02AD.f MB03QX.f select.f MB01YD.f MB01ZD.f \
- SB03OT.f MB04OY.f MB03QD.f MB04ND.f MB03QY.f \
- SB03OR.f SB03OY.f SB04PX.f MB04NY.f SB03OV.f
-
-mkoctfile SB16BD.f AB09AD.f AB09BD.f SB08GD.f SB08HD.f \
- TB01ID.f AB09AX.f MA02GD.f AB09BX.f TB01WD.f \
- MA02DD.f MB03UD.f select.f AB09DD.f SB03OU.f \
- MA02AD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \
- SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f
-
-mkoctfile SB16CD.f SB16CY.f AB09IX.f SB03OD.f MB02UD.f \
- AB09DD.f MA02AD.f MB03UD.f select.f SB03OU.f \
- MB01SD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \
- SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f
-
-system ("rm *.o");
-cd (homedir);
-
Deleted: trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m 2011-10-15 19:22:25 UTC (rev 8751)
+++ trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m 2011-10-16 09:47:01 UTC (rev 8752)
@@ -1,35 +0,0 @@
-homedir = pwd ();
-develdir = fileparts (which ("makefile_slmodred"));
-srcdir = [develdir, "/../src"];
-cd (srcdir);
-
-mkoctfile AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \
- AB09IX.f AB09IY.f SB08CD.f MB04ND.f TB01XD.f \
- MB04OD.f MB01WD.f MB03UD.f AB07MD.f SB01FY.f \
- AB09DD.f TB01LD.f SB03OU.f TB01UD.f MA02AD.f \
- MA02BD.f MB03OY.f MB03QX.f MB01PD.f select.f \
- MB01YD.f MB04NY.f MB01ZD.f SB03OT.f MB04OX.f \
- MB04OY.f MB03QD.f SB03OY.f MB03QY.f MB01QD.f \
- SB03OR.f SB03OV.f SB04PX.f
-
-mkoctfile "-Wl,-framework" "-Wl,vecLib" \
- slab09jd.cc \
- AB09JD.f TB01ID.f TB01KD.f AB07ND.f AB09JV.f \
- AB09JW.f AB09CX.f AG07BD.f AB08MD.f AB04MD.f \
- TB01LD.f delctg.f SB04PY.f AB09AX.f AB08NX.f \
- MB01SD.f AB09JX.f MA02AD.f TB01WD.f MB03OY.f \
- MB03PY.f MA02DD.f MB03UD.f MB03QX.f select.f \
- SB04PX.f SB03OU.f MB03QD.f MB03QY.f SB03OT.f \
- MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \
- MB04OY.f SB03OV.f
-
-mkoctfile AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \
- AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \
- SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \
- SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \
- SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \
- SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f
-
-system ("rm *.o");
-cd (homedir);
-
Added: trunk/octave-forge/extra/control-devel/src/IB01AD.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/IB01AD.f (rev 0)
+++ trunk/octave-forge/extra/control-devel/src/IB01AD.f 2011-10-16 09:47:01 UTC (rev 8752)
@@ -0,0 +1,686 @@
+ SUBROUTINE IB01AD( METH, ALG, JOBD, BATCH, CONCT, CTRL, NOBR, M,
+ $ L, NSMP, U, LDU, Y, LDY, N, R, LDR, SV, RCOND,
+ $ TOL, IWORK, DWORK, LDWORK, IWARN, INFO )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To preprocess the input-output data for estimating the matrices
+C of a linear time-invariant dynamical system and to find an
+C estimate of the system order. The input-output data can,
+C optionally, be processed sequentially.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C METH CHARACTER*1
+C Specifies the subspace identification method to be used,
+C as follows:
+C = 'M': MOESP algorithm with past inputs and outputs;
+C = 'N': N4SID algorithm.
+C
+C ALG CHARACTER*1
+C Specifies the algorithm for computing the triangular
+C factor R, as follows:
+C = 'C': Cholesky algorithm applied to the correlation
+C matrix of the input-output data;
+C = 'F': Fast QR algorithm;
+C = 'Q': QR algorithm applied to the concatenated block
+C Hankel matrices.
+C
+C JOBD CHARACTER*1
+C Specifies whether or not the matrices B and D should later
+C be computed using the MOESP approach, as follows:
+C = 'M': the matrices B and D should later be computed
+C using the MOESP approach;
+C = 'N': the matrices B and D should not be computed using
+C the MOESP approach.
+C This parameter is not relevant for METH = 'N'.
+C
+C BATCH CHARACTER*1
+C Specifies whether or not sequential data processing is to
+C be used, and, for sequential processing, whether or not
+C the current data block is the first block, an intermediate
+C block, or the last block, as follows:
+C = 'F': the first block in sequential data processing;
+C = 'I': an intermediate block in sequential data
+C processing;
+C = 'L': the last block in sequential data processing;
+C = 'O': one block only (non-sequential data processing).
+C NOTE that when 100 cycles of sequential data processing
+C are completed for BATCH = 'I', a warning is
+C issued, to prevent for an infinite loop.
+C
+C CONCT CHARACTER*1
+C Specifies whether or not the successive data blocks in
+C sequential data processing belong to a single experiment,
+C as follows:
+C = 'C': the current data block is a continuation of the
+C previous data block and/or it will be continued
+C by the next data block;
+C = 'N': there is no connection between the current data
+C block and the previous and/or the next ones.
+C This parameter is not used if BATCH = 'O'.
+C
+C CTRL CHARACTER*1
+C Specifies whether or not the user's confirmation of the
+C system order estimate is desired, as follows:
+C = 'C': user's confirmation;
+C = 'N': no confirmation.
+C If CTRL = 'C', a reverse communication routine, IB01OY,
+C is indirectly called (by SLICOT Library routine IB01OD),
+C and, after inspecting the singular values and system order
+C estimate, n, the user may accept n or set a new value.
+C IB01OY is not called if CTRL = 'N'.
+C
+C Input/Output Parameters
+C
+C NOBR (input) INTEGER
+C The number of block rows, s, in the input and output
+C block Hankel matrices to be processed. NOBR > 0.
+C (In the MOESP theory, NOBR should be larger than n,
+C the estimated dimension of state vector.)
+C
+C M (input) INTEGER
+C The number of system inputs. M >= 0.
+C When M = 0, no system inputs are processed.
+C
+C L (input) INTEGER
+C The number of system outputs. L > 0.
+C
+C NSMP (input) INTEGER
+C The number of rows of matrices U and Y (number of
+C samples, t). (When sequential data processing is used,
+C NSMP is the number of samples of the current data
+C block.)
+C NSMP >= 2*(M+L+1)*NOBR - 1, for non-sequential
+C processing;
+C NSMP >= 2*NOBR, for sequential processing.
+C The total number of samples when calling the routine with
+C BATCH = 'L' should be at least 2*(M+L+1)*NOBR - 1.
+C The NSMP argument may vary from a cycle to another in
+C sequential data processing, but NOBR, M, and L should
+C be kept constant. For efficiency, it is advisable to use
+C NSMP as large as possible.
+C
+C U (input) DOUBLE PRECISION array, dimension (LDU,M)
+C The leading NSMP-by-M part of this array must contain the
+C t-by-m input-data sequence matrix U,
+C U = [u_1 u_2 ... u_m]. Column j of U contains the
+C NSMP values of the j-th input component for consecutive
+C time increments.
+C If M = 0, this array is not referenced.
+C
+C LDU INTEGER
+C The leading dimension of the array U.
+C LDU >= NSMP, if M > 0;
+C LDU >= 1, if M = 0.
+C
+C Y (input) DOUBLE PRECISION array, dimension (LDY,L)
+C The leading NSMP-by-L part of this array must contain the
+C t-by-l output-data sequence matrix Y,
+C Y = [y_1 y_2 ... y_l]. Column j of Y contains the
+C NSMP values of the j-th output component for consecutive
+C time increments.
+C
+C LDY INTEGER
+C The leading dimension of the array Y. LDY >= NSMP.
+C
+C N (output) INTEGER
+C The estimated order of the system.
+C If CTRL = 'C', the estimated order has been reset to a
+C value specified by the user.
+C
+C R (output or input/output) DOUBLE PRECISION array, dimension
+C ( LDR,2*(M+L)*NOBR )
+C On exit, if ALG = 'C' and BATCH = 'F' or 'I', the leading
+C 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular part of this
+C array contains the current upper triangular part of the
+C correlation matrix in sequential data processing.
+C If ALG = 'F' and BATCH = 'F' or 'I', the array R is not
+C referenced.
+C On exit, if INFO = 0, ALG = 'Q', and BATCH = 'F' or 'I',
+C the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular
+C part of this array contains the current upper triangular
+C factor R from the QR factorization of the concatenated
+C block Hankel matrices. Denote R_ij, i,j = 1:4, the
+C ij submatrix of R, partitioned by M*NOBR, M*NOBR,
+C L*NOBR, and L*NOBR rows and columns.
+C On exit, if INFO = 0 and BATCH = 'L' or 'O', the leading
+C 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular part of
+C this array contains the matrix S, the processed upper
+C triangular factor R from the QR factorization of the
+C concatenated block Hankel matrices, as required by other
+C subroutines. Specifically, let S_ij, i,j = 1:4, be the
+C ij submatrix of S, partitioned by M*NOBR, L*NOBR,
+C M*NOBR, and L*NOBR rows and columns. The submatrix
+C S_22 contains the matrix of left singular vectors needed
+C subsequently. Useful information is stored in S_11 and
+C in the block-column S_14 : S_44. For METH = 'M' and
+C JOBD = 'M', the upper triangular part of S_31 contains
+C the upper triangular factor in the QR factorization of the
+C matrix R_1c = [ R_12' R_22' R_11' ]', and S_12
+C contains the corresponding leading part of the transformed
+C matrix R_2c = [ R_13' R_23' R_14' ]'. For METH = 'N',
+C the subarray S_41 : S_43 contains the transpose of the
+C matrix contained in S_14 : S_34.
+C The details of the contents of R need not be known if this
+C routine is followed by SLICOT Library routine IB01BD.
+C On entry, if ALG = 'C', or ALG = 'Q', and BATCH = 'I' or
+C 'L', the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper
+C triangular part of this array must contain the upper
+C triangular matrix R computed at the previous call of this
+C routine in sequential data processing. The array R need
+C not be set on entry if ALG = 'F' or if BATCH = 'F' or 'O'.
+C
+C LDR INTEGER
+C The leading dimension of the array R.
+C LDR >= MAX( 2*(M+L)*NOBR, 3*M*NOBR ),
+C for METH = 'M' and JOBD = 'M';
+C LDR >= 2*(M+L)*NOBR, for METH = 'M' and JOBD = 'N' or
+C for METH = 'N'.
+C
+C SV (output) DOUBLE PRECISION array, dimension ( L*NOBR )
+C The singular values used to estimate the system order.
+C
+C Tolerances
+C
+C RCOND DOUBLE PRECISION
+C The tolerance to be used for estimating the rank of
+C matrices. If the user sets RCOND > 0, the given value
+C of RCOND is used as a lower bound for the reciprocal
+C condition number; an m-by-n matrix whose estimated
+C condition number is less than 1/RCOND is considered to
+C be of full rank. If the user sets RCOND <= 0, then an
+C implicitly computed, default tolerance, defined by
+C RCONDEF = m*n*EPS, is used instead, where EPS is the
+C relative machine precision (see LAPACK Library routine
+C DLAMCH).
+C This parameter is not used for METH = 'M'.
+C
+C TOL DOUBLE PRECISION
+C Absolute tolerance used for determining an estimate of
+C the system order. If TOL >= 0, the estimate is
+C indicated by the index of the last singular value greater
+C than or equal to TOL. (Singular values less than TOL
+C are considered as zero.) When TOL = 0, an internally
+C computed default value, TOL = NOBR*EPS*SV(1), is used,
+C where SV(1) is the maximal singular value, and EPS is
+C the relative machine precision (see LAPACK Library routine
+C DLAMCH). When TOL < 0, the estimate is indicated by the
+C index of the singular value that has the largest
+C logarithmic gap to its successor.
+C
+C Workspace
+C
+C IWORK INTEGER array, dimension (LIWORK)
+C LIWORK >= (M+L)*NOBR, if METH = 'N';
+C LIWORK >= M+L, if METH = 'M' and ALG = 'F';
+C LIWORK >= 0, if METH = 'M' and ALG = 'C' or 'Q'.
+C
+C DWORK DOUBLE PRECISION array, dimension (LDWORK)
+C On exit, if INFO = 0, DWORK(1) returns the optimal value
+C of LDWORK, and, for METH = 'N', and BATCH = 'L' or
+C 'O', DWORK(2) and DWORK(3) contain the reciprocal
+C condition numbers of the triangular factors of the
+C matrices U_f and r_1 [6].
+C On exit, if INFO = -23, DWORK(1) returns the minimum
+C value of LDWORK.
+C Let
+C k = 0, if CONCT = 'N' and ALG = 'C' or 'Q';
+C k = 2*NOBR-1, if CONCT = 'C' and ALG = 'C' or 'Q';
+C k = 2*NOBR*(M+L+1), if CONCT = 'N' and ALG = 'F';
+C k = 2*NOBR*(M+L+2), if CONCT = 'C' and ALG = 'F'.
+C The first (M+L)*k elements of DWORK should be preserved
+C during successive calls of the routine with BATCH = 'F'
+C or 'I', till the final call with BATCH = 'L'.
+C
+C LDWORK INTEGER
+C The length of the array DWORK.
+C LDWORK >= (4*NOBR-2)*(M+L), if ALG = 'C', BATCH = 'F' or
+C 'I' and CONCT = 'C';
+C LDWORK >= 1, if ALG = 'C', BATCH = 'F' or 'I' and
+C CONCT = 'N';
+C LDWORK >= max((4*NOBR-2)*(M+L), 5*L*NOBR), if METH = 'M',
+C ALG = 'C', BATCH = 'L' and CONCT = 'C';
+C LDWORK >= max((2*M-1)*NOBR, (M+L)*NOBR, 5*L*NOBR),
+C if METH = 'M', JOBD = 'M', ALG = 'C',
+C BATCH = 'O', or
+C (BATCH = 'L' and CONCT = 'N');
+C LDWORK >= 5*L*NOBR, if METH = 'M', JOBD = 'N', ALG = 'C',
+C BATCH = 'O', or
+C (BATCH = 'L' and CONCT = 'N');
+C LDWORK >= 5*(M+L)*NOBR+1, if METH = 'N', ALG = 'C', and
+C BATCH = 'L' or 'O';
+C LDWORK >= (M+L)*2*NOBR*(M+L+3), if ALG = 'F',
+C BATCH <> 'O' and CONCT = 'C';
+C LDWORK >= (M+L)*2*NOBR*(M+L+1), if ALG = 'F',
+C BATCH = 'F', 'I' and CONCT = 'N';
+C LDWORK >= (M+L)*4*NOBR*(M+L+1)+(M+L)*2*NOBR, if ALG = 'F',
+C BATCH = 'L' and CONCT = 'N', or
+C BATCH = 'O';
+C LDWORK >= 4*(M+L)*NOBR, if ALG = 'Q', BATCH = 'F', and
+C LDR >= NS = NSMP - 2*NOBR + 1;
+C LDWORK >= max(4*(M+L)*NOBR, 5*L*NOBR), if METH = 'M',
+C ALG = 'Q', BATCH = 'O', and LDR >= NS;
+C LDWORK >= 5*(M+L)*NOBR+1, if METH = 'N', ALG = 'Q',
+C BATCH = 'O', and LDR >= NS;
+C LDWORK >= 6*(M+L)*NOBR, if ALG = 'Q', (BATCH = 'F' or 'O',
+C and LDR < NS), or (BATCH = 'I' or
+C 'L' and CONCT = 'N');
+C LDWORK >= 4*(NOBR+1)*(M+L)*NOBR, if ALG = 'Q', BATCH = 'I'
+C or 'L' and CONCT = 'C'.
+C The workspace used for ALG = 'Q' is
+C LDRWRK*2*(M+L)*NOBR + 4*(M+L)*NOBR,
+C where LDRWRK = LDWORK/(2*(M+L)*NOBR) - 2; recommended
+C value LDRWRK = NS, assuming a large enough cache size.
+C For good performance, LDWORK should be larger.
+C
+C Warning Indicator
+C
+C IWARN INTEGER
+C = 0: no warning;
+C = 1: the number of 100 cycles in sequential data
+C processing has been exhausted without signaling
+C that the last block of data was get; the cycle
+C counter was reinitialized;
+C = 2: a fast algorithm was requested (ALG = 'C' or 'F'),
+C but it failed, and the QR algorithm was then used
+C (non-sequential data processing);
+C = 3: all singular values were exactly zero, hence N = 0
+C (both input and output were identically zero);
+C = 4: the least squares problems with coefficient matrix
+C U_f, used for computing the weighted oblique
+C projection (for METH = 'N'), have a rank-deficient
+C coefficient matrix;
+C = 5: the least squares problem with coefficient matrix
+C r_1 [6], used for computing the weighted oblique
+C projection (for METH = 'N'), has a rank-deficient
+C coefficient matrix.
+C NOTE: the values 4 and 5 of IWARN have no significance
+C for the identification problem.
+C
+C Error Indicator
+C
+C INFO INTEGER
+C = 0: successful exit;
+C < 0: if INFO = -i, the i-th argument had an illegal
+C value;
+C = 1: a fast algorithm was requested (ALG = 'C', or 'F')
+C in sequential data processing, but it failed; the
+C routine can be repeatedly called again using the
+C standard QR algorithm;
+C = 2: the singular value decomposition (SVD) algorithm did
+C not converge.
+C
+C METHOD
+C
+C The procedure consists in three main steps, the first step being
+C performed by one of the three algorithms included.
+C
+C 1.a) For non-sequential data processing using QR algorithm, a
+C t x 2(m+l)s matrix H is constructed, where
+C
+C H = [ Uf' Up' Y' ], for METH = 'M',
+C s+1,2s,t 1,s,t 1,2s,t
+C
+C H = [ U' Y' ], for METH = 'N',
+C 1,2s,t 1,2s,t
+C
+C and Up , Uf , U , and Y are block Hankel
+C 1,s,t s+1,2s,t 1,2s,t 1,2s,t
+C matrices defined in terms of the input and output data [3].
+C A QR factorization is used to compress the data.
+C The fast QR algorithm uses a QR factorization which exploits
+C the block-Hankel structure. Actually, the Cholesky factor of H'*H
+C is computed.
+C
+C 1.b) For sequential data processing using QR algorithm, the QR
+C decomposition is done sequentially, by updating the upper
+C triangular factor R. This is also performed internally if the
+C workspace is not large enough to accommodate an entire batch.
+C
+C 1.c) For non-sequential or sequential data processing using
+C Cholesky algorithm, the correlation matrix of input-output data is
+C computed (sequentially, if requested), taking advantage of the
+C block Hankel structure [7]. Then, the Cholesky factor of the
+C correlation matrix is found, if possible.
+C
+C 2) A singular value decomposition (SVD) of a certain matrix is
+C then computed, which reveals the order n of the system as the
+C number of "non-zero" singular values. For the MOESP approach, this
+C matrix is [ R_24' R_34' ]' := R(ms+1:(2m+l)s,(2m+l)s+1:2(m+l)s),
+C where R is the upper triangular factor R constructed by SLICOT
+C Library routine IB01MD. For the N4SID approach, a weighted
+C oblique projection is computed from the upper triangular factor R
+C and its SVD is then found.
+C
+C 3) The singular values are compared to the given, or default TOL,
+C and the estimated order n is returned, possibly after user's
+C confirmation.
+C
+C REFERENCES
+C
+C [1] Verhaegen M., and Dewilde, P.
+C Subspace Model Identification. Part 1: The output-error
+C state-space model identification class of algorithms.
+C Int. J. Control, 56, pp. 1187-1210, 1992.
+C
+C [2] Verhaegen M.
+C Subspace Model Identification. Part 3: Analysis of the
+C ordinary output-error state-space model identification
+C algorithm.
+C Int. J. Control, 58, pp. 555-586, 1993.
+C
+C [3] Verhaegen M.
+C Identification of the deterministic part of MIMO state space
+C models given in innovations form from input-output data.
+C Automatica, Vol.30, No.1, pp.61-74, 1994.
+C
+C [4] Van Overschee, P., and De Moor, B.
+C N4SID: Subspace Algorithms for the Identification of
+C Combined Deterministic-Stochastic Systems.
+C Automatica, Vol.30, No.1, pp. 75-93, 1994.
+C
+C [5] Peternell, K., Scherrer, W. and Deistler, M.
+C Statistical Analysis of Novel Subspace Identification Methods.
+C Signal Processing, 52, pp. 161-177, 1996.
+C
+C [6] Sima, V.
+C Subspace-based Algorithms for Multivariable System
+C Identification.
+C Studies in Informatics and Control, 5, pp. 335-344, 1996.
+C
+C [7] Sima, V.
+C Cholesky or QR Factorization for Data Compression in
+C Subspace-based Identification ?
+C Proceedings of the Second NICONET Workshop on ``Numerical
+C Control Software: SLICOT, a Useful Tool in Industry'',
+C December 3, 1999, INRIA Rocquencourt, France, pp. 75-80, 1999.
+C
+C NUMERICAL ASPECTS
+C
+C The implemented method is numerically stable (when QR algorithm is
+C used), reliable and efficient. The fast Cholesky or QR algorithms
+C are more efficient, but the accuracy could diminish by forming the
+C correlation matrix.
+C The most time-consuming computational step is step 1:
+C 2
+C The QR algorithm needs 0(t(2(m+l)s) ) floating point operations.
+C 2 3
+C The Cholesky algorithm needs 0(2t(m+l) s)+0((2(m+l)s) ) floating
+C point operations.
+C 2 3 2
+C The fast QR algorithm needs 0(2t(m+l) s)+0(4(m+l) s ) floating
+C point operations.
+C 3
+C Step 2 of the algorithm requires 0(((m+l)s) ) floating point
+C operations.
+C
+C FURTHER COMMENTS
+C
+C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the
+C calculations could be rather inefficient if only minimal workspace
+C (see argument LDWORK) is provided. It is advisable to provide as
+C much workspace as possible. Almost optimal efficiency can be
+C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the
+C cache size is large enough to accommodate R, U, Y, and DWORK.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Katholieke Universiteit Leuven, Feb. 2000.
+C
+C REVISIONS
+C
+C August 2000, March 2005.
+C
+C KEYWORDS
+C
+C Cholesky decomposition, Hankel matrix, identification methods,
+C multivariable systems, QR decomposition, singular value
+C decomposition.
+C
+C ******************************************************************
+C
+C .. Scalar Arguments ..
+ DOUBLE PRECISION RCOND, TOL
+ INTEGER INFO, IWARN, L, LDR, LDU, LDWORK, LDY, M, N,
+ $ NOBR, NSMP
+ CHARACTER ALG, BATCH, CONCT, CTRL, JOBD, METH
+C .. Array Arguments ..
+ INTEGER IWORK(*)
+ DOUBLE PRECISION DWORK(*), R(LDR, *), SV(*), U(LDU, *),
+ $ Y(LDY, *)
+C .. Local Scalars ..
+ INTEGER IWARNL, LMNOBR, LNOBR, MAXWRK, MINWRK, MNOBR,
+ $ NOBR21, NR, NS, NSMPSM
+ LOGICAL CHALG, CONNEC, CONTRL, FIRST, FQRALG, INTERM,
+ $ JOBDM, LAST, MOESP, N4SID, ONEBCH, QRALG
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. External Subroutines ..
+ EXTERNAL IB01MD, IB01ND, IB01OD, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC MAX
+C .. Save Statement ..
+C MAXWRK is used to store the optimal workspace.
+C NSMPSM is used to sum up the NSMP values for BATCH <> 'O'.
+ SAVE MAXWRK, NSMPSM
+C ..
+C .. Executable Statements ..
+C
+C Decode the scalar input parameters.
+C
+ MOESP = LSAME( METH, 'M' )
+ N4SID = LSAME( METH, 'N' )
+ FQRALG = LSAME( ALG, 'F' )
+ QRALG = LSAME( ALG, 'Q' )
+ CHALG = LSAME( ALG, 'C' )
+ JOBDM = LSAME( JOBD, 'M' )
+ ONEBCH = LSAME( BATCH, 'O' )
+ FIRST = LSAME( BATCH, 'F' ) .OR. ONEBCH
+ INTERM = LSAME( BATCH, 'I' )
+ LAST = LSAME( BATCH, 'L' ) .OR. ONEBCH
+ CONTRL = LSAME( CTRL, 'C' )
+C
+ IF( .NOT.ONEBCH ) THEN
+ CONNEC = LSAME( CONCT, 'C' )
+ ELSE
+ CONNEC = .FALSE.
+ END IF
+C
+ MNOBR = M*NOBR
+ LNOBR = L*NOBR
+ LMNOBR = LNOBR + MNOBR
+ NR = LMNOBR + LMNOBR
+ NOBR21 = 2*NOBR - 1
+ IWARN = 0
+ INFO = 0
+ IF( FIRST ) THEN
+ MAXWRK = 1
+ NSMPSM = 0
+ END IF
+ NSMPSM = NSMPSM + NSMP
+C
+C Check the scalar input parameters.
+C
+ IF( .NOT.( MOESP .OR. N4SID ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.( FQRALG .OR. QRALG .OR. CHALG ) ) THEN
+ INFO = -2
+ ELSE IF( MOESP .AND. .NOT.( JOBDM .OR. LSAME( JOBD, 'N' ) ) ) THEN
+ INFO = -3
+ ELSE IF( .NOT.( FIRST .OR. INTERM .OR. LAST ) ) THEN
+ INFO = -4
+ ELSE IF( .NOT. ONEBCH ) THEN
+ IF( .NOT.( CONNEC .OR. LSAME( CONCT, 'N' ) ) )
+ $ INFO = -5
+ END IF
+ IF( INFO.EQ.0 ) THEN
+ IF( .NOT.( CONTRL .OR. LSAME( CTRL, 'N' ) ) ) THEN
+ INFO = -6
+ ELSE IF( NOBR.LE.0 ) THEN
+ INFO = -7
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -8
+ ELSE IF( L.LE.0 ) THEN
+ INFO = -9
+ ELSE IF( NSMP.LT.2*NOBR .OR.
+ $ ( LAST .AND. NSMPSM.LT.NR+NOBR21 ) ) THEN
+ INFO = -10
+ ELSE IF( LDU.LT.1 .OR. ( M.GT.0 .AND. LDU.LT.NSMP ) ) THEN
+ INFO = -12
+ ELSE IF( LDY.LT.NSMP ) THEN
+ INFO = -14
+ ELSE IF( LDR.LT.NR .OR. ( MOESP .AND. JOBDM .AND.
+ $ LDR.LT.3*MNOBR ) ) THEN
+ INFO = -17
+ ELSE
+C
+C Compute workspace.
+C (Note: Comments in the code beginning "Workspace:" describe
+C the minimal amount of workspace needed at that point in the
+C code, as well as the preferred amount for good performance.)
+C
+ NS = NSMP - NOBR21
+ IF ( CHALG ) THEN
+ IF ( .NOT.LAST ) THEN
+ IF ( CONNEC ) THEN
+ MINWRK = 2*( NR - M - L )
+ ELSE
+ MINWRK = 1
+ END IF
+ ELSE IF ( MOESP ) THEN
+ IF ( CONNEC .AND. .NOT.ONEBCH ) THEN
+ MINWRK = MAX( 2*( NR - M - L ), 5*LNOBR )
+ ELSE
+ MINWRK = 5*LNOBR
+ IF ( JOBDM )
+ $ MINWRK = MAX( 2*MNOBR - NOBR, LMNOBR, MINWRK )
+ END IF
+ ELSE
+ MINWRK = 5*LMNOBR + 1
+ END IF
+ ELSE IF ( FQRALG ) THEN
+ IF ( .NOT.ONEBCH .AND. CONNEC ) THEN
+ MINWRK = NR*( M + L + 3 )
+ ELSE IF ( FIRST .OR. INTERM ) THEN
+ MINWRK = NR*( M + L + 1 )
+ ELSE
+ MINWRK = 2*NR*( M + L + 1 ) + NR
+ END IF
+ ELSE
+ MINWRK = 2*NR
+ IF ( ONEBCH .AND. LDR.GE.NS ) THEN
+ IF ( MOESP ) THEN
+ MINWRK = MAX( MINWRK, 5*LNOBR )
+ ELSE
+ MINWRK = 5*LMNOBR + 1
+ END IF
+ END IF
+ IF ( FIRST ) THEN
+ IF ( LDR.LT.NS ) THEN
+ MINWRK = MINWRK + NR
+ END IF
+ ELSE
+ IF ( CONNEC ) THEN
+ MINWRK = MINWRK*( NOBR + 1 )
+ ELSE
+ MINWRK = MINWRK + NR
+ END IF
+ END IF
+ END IF
+C
+ MAXWRK = MINWRK
+C
+ IF( LDWORK.LT.MINWRK ) THEN
+ INFO = -23
+ DWORK( 1 ) = MINWRK
+ END IF
+ END IF
+ END IF
+C
+C Return if there are illegal arguments.
+C
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'IB01AD', -INFO )
+ RETURN
+ END IF
+C
+C Compress the input-output data.
+C Workspace: need c*(M+L)*NOBR, where c is a constant depending
+C on the algorithm and the options used
+C (see SLICOT Library routine IB01MD);
+C prefer larger.
+C
+ CALL IB01MD( METH, ALG, BATCH, CONCT, NOBR, M, L, NSMP, U, LDU, Y,
+ $ LDY, R, LDR, IWORK, DWORK, LDWORK, IWARN, INFO )
+C
+ IF ( INFO.EQ.1 ) THEN
+C
+C Error return: A fast algorithm was requested (ALG = 'C', 'F')
+C in sequential data processing, but it failed.
+C
+ RETURN
+ END IF
+C
+ MAXWRK = MAX( MAXWRK, INT( DWORK( 1 ) ) )
+C
+ IF ( .NOT.LAST ) THEN
+C
+C Return to get new data.
+C
+ RETURN
+ END IF
+C
+C Find the singular value decomposition (SVD) giving the system
+C order, and perform related preliminary calculations needed for
+C computing the system matrices.
+C Workspace: need max( (2*M-1)*NOBR, (M+L)*NOBR, 5*L*NOBR ),
+C if METH = 'M';
+C 5*(M+L)*NOBR+1, if METH = 'N';
+C prefer larger.
+C
+ CALL IB01ND( METH, JOBD, NOBR, M, L, R, LDR, SV, RCOND, IWORK,
+ $ DWORK, LDWORK, IWARNL, INFO )
+ IWARN = MAX( IWARN, IWARNL )
+C
+ IF ( INFO.EQ.2 ) THEN
+C
+C Error return: the singular value decomposition (SVD) algorithm
+C did not converge.
+C
+ RETURN
+ END IF
+C
+C Estimate the system order.
+C
+ CALL IB01OD( CTRL, NOBR, L, SV, N, TOL, IWARNL, INFO )
+ IWARN = MAX( IWARN, IWARNL )
+C
+C Return optimal workspace in DWORK(1).
+C
+ DWORK( 1 ) = MAX( MAXWRK, INT( DWORK( 1 ) ) )
+ RETURN
+C
+C *** Last line of IB01AD ***
+ END
Property changes on: trunk/octave-forge/extra/control-devel/src/IB01AD.f
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/control-devel/src/IB01BD.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/IB01BD.f (rev 0)
+++ trunk/octave-forge/extra/control-devel/src/IB01BD.f 2011-10-16 09:47:01 UTC (rev 8752)
@@ -0,0 +1,791 @@
+ SUBROUTINE IB01BD( METH, JOB, JOBCK, NOBR, N, M, L, NSMPL, R,
+ $ LDR, A, LDA, C, LDC, B, LDB, D, LDD, Q, LDQ,
+ $ RY, LDRY, S, LDS, K, LDK, TOL, IWORK, DWORK,
+ $ LDWORK, BWORK, IWARN, INFO )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To estimate the system matrices A, C, B, and D, the noise
+C covariance matrices Q, Ry, and S, and the Kalman gain matrix K
+C of a linear time-invariant state space model, using the
+C processed triangular factor R of the concatenated block Hankel
+C matrices, provided by SLICOT Library routine IB01AD.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C METH CHARACTER*1
+C Specifies the subspace identification method to be used,
+C as follows:
+C = 'M': MOESP algorithm with past inputs and outputs;
+C = 'N': N4SID algorithm;
+C = 'C': combined method: MOESP algorithm for finding the
+C matrices A and C, and N4SID algorithm for
+C finding the matrices B and D.
+C
+C JOB CHARACTER*1
+C Specifies which matrices should be computed, as follows:
+C = 'A': compute all system matrices, A, B, C, and D;
+C = 'C': compute the matrices A and C only;
+C = 'B': compute the matrix B only;
+C = 'D': compute the matrices B and D only.
+C
+C JOBCK CHARACTER*1
+C Specifies whether or not the covariance matrices and the
+C Kalman gain matrix are to be computed, as follows:
+C = 'C': the covariance matrices only should be computed;
+C = 'K': the covariance matrices and the Kalman gain
+C matrix should be computed;
+C = 'N': the covariance matrices and the Kalman gain matrix
+C should not be computed.
+C
+C Input/Output Parameters
+C
+C NOBR (input) INTEGER
+C The number of block rows, s, in the input and output
+C Hankel matrices processed by other routines. NOBR > 1.
+C
+C N (input) INTEGER
+C The order of the system. NOBR > N > 0.
+C
+C M (input) INTEGER
+C The number of system inputs. M >= 0.
+C
+C L (input) INTEGER
+C The number of system outputs. L > 0.
+C
+C NSMPL (input) INTEGER
+C If JOBCK = 'C' or 'K', the total number of samples used
+C for calculating the covariance matrices.
+C NSMPL >= 2*(M+L)*NOBR.
+C This parameter is not meaningful if JOBCK = 'N'.
+C
+C R (input/workspace) DOUBLE PRECISION array, dimension
+C ( LDR,2*(M+L)*NOBR )
+C On entry, the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR part
+C of this array must contain the relevant data for the MOESP
+C or N4SID algorithms, as constructed by SLICOT Library
+C routine IB01AD. Let R_ij, i,j = 1:4, be the
+C ij submatrix of R (denoted S in IB01AD), partitioned
+C by M*NOBR, L*NOBR, M*NOBR, and L*NOBR rows and
+C columns. The submatrix R_22 contains the matrix of left
+C singular vectors used. Also needed, for METH = 'N' or
+C JOBCK <> 'N', are the submatrices R_11, R_14 : R_44,
+C and, for METH = 'M' or 'C' and JOB <> 'C', the
+C submatrices R_31 and R_12, containing the processed
+C matrices R_1c and R_2c, respectively, as returned by
+C SLICOT Library routine IB01AD.
+C Moreover, if METH = 'N' and JOB = 'A' or 'C', the
+C block-row R_41 : R_43 must contain the transpose of the
+C block-column R_14 : R_34 as returned by SLICOT Library
+C routine IB01AD.
+C The remaining part of R is used as workspace.
+C On exit, part of this array is overwritten. Specifically,
+C if METH = 'M', R_22 and R_31 are overwritten if
+C JOB = 'B' or 'D', and R_12, R_22, R_14 : R_34,
+C and possibly R_11 are overwritten if JOBCK <> 'N';
+C if METH = 'N', all needed submatrices are overwritten.
+C The details of the contents of R need not be known if
+C this routine is called once just after calling the SLICOT
+C Library routine IB01AD.
+C
+C LDR INTEGER
+C The leading dimension of the array R.
+C LDR >= 2*(M+L)*NOBR.
+C
+C A (input or output) DOUBLE PRECISION array, dimension
+C (LDA,N)
+C On entry, if METH = 'N' or 'C' and JOB = 'B' or 'D',
+C the leading N-by-N part of this array must contain the
+C system state matrix.
+C If METH = 'M' or (METH = 'N' or 'C' and JOB = 'A'
+C or 'C'), this array need not be set on input.
+C On exit, if JOB = 'A' or 'C' and INFO = 0, the
+C leading N-by-N part of this array contains the system
+C state matrix.
+C
+C LDA INTEGER
+C The leading dimension of the array A.
+C LDA >= N, if JOB = 'A' or 'C', or METH = 'N' or 'C'
+C and JOB = 'B' or 'D';
+C LDA >= 1, otherwise.
+C
+C C (input or output) DOUBLE PRECISION array, dimension
+C (LDC,N)
+C On entry, if METH = 'N' or 'C' and JOB = 'B' or 'D',
+C the leading L-by-N part of this array must contain the
+C system output matrix.
+C If METH = 'M' or (METH = 'N' or 'C' and JOB = 'A'
+C or 'C'), this array need not be set on input.
+C On exit, if JOB = 'A' or 'C' and INFO = 0, or
+C INFO = 3 (or INFO >= 0, for METH = 'M'), the leading
+C L-by-N part of this array contains the system output
+C matrix.
+C
+C LDC INTEGER
+C The leading dimension of the array C.
+C LDC >= L, if JOB = 'A' or 'C', or METH = 'N' or 'C'
+C and JOB = 'B' or 'D';
+C LDC >= 1, otherwise.
+C
+C B (output) DOUBLE PRECISION array, dimension (LDB,M)
+C If M > 0, JOB = 'A', 'B', or 'D' and INFO = 0, the
+C leading N-by-M part of this array contains the system
+C input matrix. If M = 0 or JOB = 'C', this array is
+C not referenced.
+C
+C LDB INTEGER
+C The leading dimension of the array B.
+C LDB >= N, if M > 0 and JOB = 'A', 'B', or 'D';
+C LDB >= 1, if M = 0 or JOB = 'C'.
+C
+C D (output) DOUBLE PRECISION array, dimension (LDD,M)
+C If M > 0, JOB = 'A' or 'D' and INFO = 0, the leading
+C L-by-M part of this array contains the system input-output
+C matrix. If M = 0 or JOB = 'C' or 'B', this array is
+C not referenced.
+C
+C LDD INTEGER
+C The leading dimension of the array D.
+C LDD >= L, if M > 0 and JOB = 'A' or 'D';
+C LDD >= 1, if M = 0 or JOB = 'C' or 'B'.
+C
+C Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
+C If JOBCK = 'C' or 'K', the leading N-by-N part of this
+C array contains the positive semidefinite state covariance
+C matrix. If JOBCK = 'K', this matrix has been used as
+C state weighting matrix for computing the Kalman gain.
+C This parameter is not referenced if JOBCK = 'N'.
+C
+C LDQ INTEGER
+C The leading dimension of the array Q.
+C LDQ >= N, if JOBCK = 'C' or 'K';
+C LDQ >= 1, if JOBCK = 'N'.
+C
+C RY (output) DOUBLE PRECISION array, dimension (LDRY,L)
+C If JOBCK = 'C' or 'K', the leading L-by-L part of this
+C array contains the positive (semi)definite output
+C covariance matrix. If JOBCK = 'K', this matrix has been
+C used as output weighting matrix for computing the Kalman
+C gain.
+C This parameter is not referenced if JOBCK = 'N'.
+C
+C LDRY INTEGER
+C The leading dimension of the array RY.
+C LDRY >= L, if JOBCK = 'C' or 'K';
+C LDRY >= 1, if JOBCK = 'N'.
+C
+C S (output) DOUBLE PRECISION array, dimension (LDS,L)
+C If JOBCK = 'C' or 'K', the leading N-by-L part of this
+C array contains the state-output cross-covariance matrix.
+C If JOBCK = 'K', this matrix has been used as state-
+C output weighting matrix for computing the Kalman gain.
+C This parameter is not referenced if JOBCK = 'N'.
+C
+C LDS INTEGER
+C The leading dimension of the array S.
+C LDS >= N, if JOBCK = 'C' or 'K';
+C LDS >= 1, if JOBCK = 'N'.
+C
+C K (output) DOUBLE PRECISION array, dimension ( LDK,L )
+C If JOBCK = 'K', the leading N-by-L part of this array
+C contains the estimated Kalman gain matrix.
+C If JOBCK = 'C' or 'N', this array is not referenced.
+C
+C LDK INTEGER
+C The leading dimension of the array K.
+C LDK >= N, if JOBCK = 'K';
+C LDK >= 1, if JOBCK = 'C' or 'N'.
+C
+C Tolerances
+C
+C TOL DOUBLE PRECISION
+C The tolerance to be used for estimating the rank of
+C matrices. If the user sets TOL > 0, then the given value
+C of TOL is used as a lower bound for the reciprocal
+C condition number; an m-by-n matrix whose estimated
+C condition number is less than 1/TOL is considered to
+C be of full rank. If the user sets TOL <= 0, then an
+C implicitly computed, default tolerance, defined by
+C TOLDEF = m*n*EPS, is used instead, where EPS is the
+C relative machine precision (see LAPACK Library routine
+C DLAMCH).
+C
+C Workspace
+C
+C IWORK INTEGER array, dimension (LIWORK)
+C LIWORK >= max(LIW1,LIW2), where
+C LIW1 = N, if METH <> 'N' and M = 0
+C or JOB = 'C' and JOBCK = 'N';
+C LIW1 = M*NOBR+N, if METH <> 'N', JOB = 'C',
+C and JOBCK <> 'N';
+C LIW1 = max(L*NOBR,M*NOBR), if METH = 'M', JOB <> 'C',
+C and JOBCK = 'N';
+C LIW1 = max(L*NOBR,M*NOBR+N), if METH = 'M', JOB <> 'C',
+C and JOBCK = 'C' or 'K';
+C LIW1 = max(M*NOBR+N,M*(N+L)), if METH = 'N', or METH = 'C'
+C and JOB <> 'C';
+C LIW2 = 0, if JOBCK <> 'K';
+C LIW2 = N*N, if JOBCK = 'K'.
+C
+C DWORK DOUBLE PRECISION array, dimension (LDWORK)
+C On exit, if INFO = 0, DWORK(1) returns the optimal value
+C of LDWORK, and DWORK(2), DWORK(3), DWORK(4), and
+C DWORK(5) contain the reciprocal condition numbers of the
+C triangular factors of the following matrices (defined in
+C SLICOT Library routine IB01PD and in the lower level
+C routines):
+C GaL (GaL = Un(1:(s-1)*L,1:n)),
+C R_1c (if METH = 'M' or 'C'),
+C M (if JOBCK = 'C' or 'K' or METH = 'N'), and
+C Q or T (see SLICOT Library routine IB01PY or IB01PX),
+C respectively.
+C If METH = 'N', DWORK(3) is set to one without any
+C calculations. Similarly, if METH = 'M' and JOBCK = 'N',
+C DWORK(4) is set to one. If M = 0 or JOB = 'C',
+C DWORK(3) and DWORK(5) are set to one.
+C If JOBCK = 'K' and INFO = 0, DWORK(6) to DWORK(13)
+C contain information about the accuracy of the results when
+C computing the Kalman gain matrix, as follows:
+C DWORK(6) - reciprocal condition number of the matrix
+C U11 of the Nth order system of algebraic
+C equations from which the solution matrix X
+C of the Riccati equation is obtained;
+C DWORK(7) - reciprocal pivot growth factor for the LU
+C factorization of the matrix U11;
+C DWORK(8) - reciprocal condition number of the matrix
+C As = A - S*inv(Ry)*C, which is inverted by
+C the standard Riccati solver;
+C DWORK(9) - reciprocal pivot growth factor for the LU
+C factorization of the matrix As;
+C DWORK(10) - reciprocal condition number of the matrix
+C Ry;
+C DWORK(11) - reciprocal condition number of the matrix
+C Ry + C*X*C';
+C DWORK(12) - reciprocal condition number for the Riccati
+C equation solution;
+C DWORK(13) - forward error bound for the Riccati
+C equation solution.
+C On exit, if INFO = -30, DWORK(1) returns the minimum
+C value of LDWORK.
+C
+C LDWORK INTEGER
+C The length of the array DWORK.
+C LDWORK >= max( LDW1,LDW2,LDW3 ), where, if METH = 'M',
+C LDW1 >= max( 2*(L*NOBR-L)*N+2*N, (L*NOBR-L)*N+N*N+7*N ),
+C if JOB = 'C' or JOB = 'A' and M = 0;
+C LDW1 >= max( 2*(L*NOBR-L)*N+N*N+7*N,
+C (L*NOBR-L)*N+N+6*M*NOBR, (L*NOBR-L)*N+N+
+C max( L+M*NOBR, L*NOBR +
+C max( 3*L*NOBR+1, M ) ) ),
+C if M > 0 and JOB = 'A', 'B', or 'D';
+C LDW2 >= 0, if JOBCK = 'N';
+C LDW2 >= L*NOBR*N+
+C max( (L*NOBR-L)*N+Aw+2*N+max(5*N,(2*M+L)*NOBR+L),
+C 4*(M*NOBR+N)+1, M*NOBR+2*N+L ),
+C if JOBCK = 'C' or 'K',
+C where Aw = N+N*N, if M = 0 or JOB = 'C';
+C Aw = 0, otherwise;
+C if METH = 'N',
+C LDW1 >= L*NOBR*N+max( (L*NOBR-L)*N+2*N+(2*M+L)*NOBR+L,
+C 2*(L*NOBR-L)*N+N*N+8*N,
+C N+4*(M*NOBR+N)+1, M*NOBR+3*N+L );
+C LDW2 >= 0, if M = 0 or JOB = 'C';
+C LDW2 >= L*NOBR*N+M*NOBR*(N+L)*(M*(N+L)+1)+
+C max( (N+L)**2, 4*M*(N+L)+1 ),
+C if M > 0 and JOB = 'A', 'B', or 'D';
+C and, if METH = 'C', LDW1 as
+C max( LDW1 for METH = 'M', JOB = 'C', LDW1 for METH = 'N'),
+C and LDW2 for METH = 'N' are used;
+C LDW3 >= 0, if JOBCK <> 'K';
+C LDW3 >= max( 4*N*N+2*N*L+L*L+max( 3*L,N*L ),
+C 14*N*N+12*N+5 ), if JOBCK = 'K'.
+C For good performance, LDWORK should be larger.
+C
+C BWORK LOGICAL array, dimension ...
[truncated message content] |
