# HG changeset patch
# User paramaniac
# Date 1331839460 0
# Node ID 0583d649c01a24f3889ead5bbaec273b3ac2a46d
# Parent 3828294802c5b6f8d50ecb5d378539652207aec9
control-devel: quicksave id draft code
diff -r 3828294802c5 -r 0583d649c01a extra/control-devel/src/slib01bd.cc
--- a/extra/control-devel/src/slib01bd.cc Thu Mar 15 17:44:49 2012 +0000
+++ b/extra/control-devel/src/slib01bd.cc Thu Mar 15 19:24:20 2012 +0000
@@ -63,7 +63,7 @@
int nargin = args.length ();
octave_value_list retval;
- if (nargin != 11)
+ if (nargin != 5)
{
print_usage ();
}
@@ -71,81 +71,37 @@
{
// arguments in
char meth;
- char alg;
- char jobd;
- char batch;
- char conct;
- char ctrl;
+ char job = 'A';
+ char jobck = 'K';
- Matrix y = args(0).matrix_value ();
- Matrix u = args(1).matrix_value ();
- int nobr = args(2).int_value ();
-
+ Matrix r = args(0).matrix_value ();
+ int nsmpl = args(1).int_value ();
+ int n = args(2).int_value ();
const int imeth = args(3).int_value ();
- const int ialg = args(4).int_value ();
- const int ijobd = args(5).int_value ();
- const int ibatch = args(6).int_value ();
- const int iconct = args(7).int_value ();
- const int ictrl = args(8).int_value ();
-
- double rcond = args(9).double_value ();
- double tol = args(10).double_value ();
-
+ double tol = args(4).double_value ();
- if (imeth == 0)
- meth = 'M';
- else
- meth = 'N';
- switch (ialg)
+ switch (imeth)
{
case 0:
- alg = 'C';
+ meth = 'M';
break;
case 1:
- alg = 'F';
+ meth = 'N';
break;
case 2:
- alg = 'Q';
+ meth = 'C';
break;
default:
- error ("slib01bd: argument 'alg' invalid");
- }
-
- if (ijobd == 0)
- jobd = 'M';
- else
- jobd = 'N';
-
- switch (ibatch)
- {
- case 0:
- batch = 'F';
- break;
- case 1:
- batch = 'I';
- break;
- case 2:
- batch = 'L';
- break;
- case 3:
- batch = 'O';
- break;
- default:
- error ("slib01bd: argument 'batch' invalid");
+ error ("slib01bd: argument 'meth' invalid");
}
- if (iconct == 0)
- conct = 'C';
- else
- conct = 'N';
-
- if (ictrl == 0)
- ctrl = 'C';
- else
- ctrl = 'N';
-
-
+ // TODO: if meth == 'C', which meth should be taken for IB01AD.f, 'M' or 'N'?
+
+ int m =
+ int nobr = r.rows () / (2*(m+l));
+
+
int m = u.columns (); // m: number of inputs
int l = y.columns (); // l: number of outputs
int nsmp = y.rows (); // nsmp: number of samples
@@ -178,57 +134,10 @@
// workspace
int liwork;
- if (meth == 'N') // if METH = 'N'
- liwork = (m+l)*nobr;
- else if (alg == 'F') // if METH = 'M' and ALG = 'F'
- liwork = m+l;
- else // if METH = 'M' and ALG = 'C' or 'Q'
- liwork = 0;
- // TODO: Handle 'k' for DWORK
int ldwork;
- ldwork = 0;
-
- if (alg == 'C' && (batch == 'F' || batch == 'I') && conct = 'C')
- ldwork = (4*nobr-2)*(m+l);
- else if (alg == 'C' && (batch == 'F' || batch == 'I') && conct = 'N')
- ldwork = 1;
- else if (meth == 'M' && alg == 'C' && batch == 'L' && conct == 'C')
- ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr);
- else if ((meth == 'M' && jobd = 'M' && alg == 'C' && batch == 'O') || (batch == 'L' && conct == 'N'))
- ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
- else if ((meth == 'M' && jobd == 'N' && alg == 'C' && batch == 'O') || (batch == 'L' && conct == 'N'))
- ldwork = 5*l*nobr;
-
- // FIXME : two times || (batch == 'L' && conct == 'N') doesn't make sense
-
-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.
diff -r 3828294802c5 -r 0583d649c01a extra/control-devel/src/slident.cc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/control-devel/src/slident.cc Thu Mar 15 19:24:20 2012 +0000
@@ -0,0 +1,466 @@
+/*
+
+Copyright (C) 2012 Lukas F. Reichlin
+
+This file is part of LTI Syncope.
+
+LTI Syncope is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+LTI Syncope is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with LTI Syncope. If not, see .
+
+SLICOT system identification
+Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V.
+
+
+Author: Lukas Reichlin
+Created: March 2012
+Version: 0.1
+
+*/
+
+#include
+#include
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ib01ad, IB01AD)
+ (char& METH, char& ALG, char& JOBD,
+ char& BATCH, char& CONCT, char& CTRL,
+ int& NOBR, int& M, int& L,
+ int& NSMP,
+ double* U, int& LDU,
+ double* Y, int& LDY,
+ int& N,
+ double* R, int& LDR,
+ double* SV,
+ double& RCOND, double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+
+ int F77_FUNC (ib01bd, IB01BD)
+ (char& METH, char& JOB, char& JOBCK,
+ int& NOBR, int& N, int& M, int& L,
+ int& NSMPL,
+ double* R, int& LDR,
+ double* A, int& LDA,
+ double* C, int& LDC,
+ double* B, int& LDB,
+ double* D, int& LDD,
+ double* Q, int& LDQ,
+ double* RY, int& LDRY,
+ double* S, int& LDS,
+ double* K, int& LDK,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ bool* BWORK,
+ int& IWARN, int& INFO);
+
+}
+
+// PKG_ADD: autoload ("slident", "devel_slicot_functions.oct");
+DEFUN_DLD (slident, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot IB01AD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 11)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char meth;
+ char alg;
+ char jobd;
+ char batch;
+ char conct;
+ char ctrl;
+
+ Matrix y = args(0).matrix_value ();
+ Matrix u = args(1).matrix_value ();
+ int nobr = args(2).int_value ();
+
+ const int imeth = args(3).int_value ();
+ const int ialg = args(4).int_value ();
+ const int ijobd = args(5).int_value ();
+ const int ibatch = args(6).int_value ();
+ const int iconct = args(7).int_value ();
+ const int ictrl = args(8).int_value ();
+
+ double rcond = args(9).double_value ();
+ double tol = args(10).double_value ();
+
+
+ if (imeth == 0)
+ meth = 'M';
+ else
+ meth = 'N';
+
+ switch (ialg)
+ {
+ case 0:
+ alg = 'C';
+ break;
+ case 1:
+ alg = 'F';
+ break;
+ case 2:
+ alg = 'Q';
+ break;
+ default:
+ error ("slib01ad: argument 'alg' invalid");
+ }
+
+ if (ijobd == 0)
+ jobd = 'M';
+ else
+ jobd = 'N';
+
+ switch (ibatch)
+ {
+ case 0:
+ batch = 'F';
+ break;
+ case 1:
+ batch = 'I';
+ break;
+ case 2:
+ batch = 'L';
+ break;
+ case 3:
+ batch = 'O';
+ break;
+ default:
+ error ("slib01ad: argument 'batch' invalid");
+ }
+
+ if (iconct == 0)
+ conct = 'C';
+ else
+ conct = 'N';
+
+ if (ictrl == 0)
+ ctrl = 'C';
+ else
+ ctrl = 'N';
+
+
+ int m = u.columns (); // m: number of inputs
+ int l = y.columns (); // l: number of outputs
+ int nsmp = y.rows (); // nsmp: number of samples
+ // y.rows == u.rows is checked by iddata class
+ // TODO: check minimal nsmp size
+
+ int ldu;
+
+ if (m == 0)
+ ldu = 1;
+ else // m > 0
+ ldu = nsmp;
+
+ int ldy = nsmp;
+
+ // arguments out
+ int n;
+ int ldr;
+
+ if (meth == 'M' && jobd == 'M')
+ ldr = max (2*(m+l)*nobr, 3*m*nobr);
+ else if (meth == 'N' || (meth == 'M' && jobd == 'N'))
+ ldr = 2*(m+l)*nobr;
+ else
+ error ("slib01ad: could not handle 'ldr' case");
+
+ Matrix r (ldr, 2*(m+l)*nobr);
+ ColumnVector sv (l*nobr);
+
+ // workspace
+ int liwork;
+
+ if (meth == 'N') // if METH = 'N'
+ liwork = (m+l)*nobr;
+ else if (alg == 'F') // if METH = 'M' and ALG = 'F'
+ liwork = m+l;
+ else // if METH = 'M' and ALG = 'C' or 'Q'
+ liwork = 0;
+
+ // TODO: Handle 'k' for DWORK
+
+ int ldwork;
+
+ ldwork = 0;
+
+ if (alg == 'C')
+ {
+ if (batch == 'F' || batch == 'I')
+ {
+ if (conct == 'C')
+ ldwork = (4*nobr-2)*(m+l);
+ else // (conct == 'N')
+ ldwork = 1;
+ }
+ else if (meth == 'M') // && (batch == 'L' || batch == 'O')
+ {
+ if (conct == 'C' && batch == 'L')
+ ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr);
+ else if (jobd == 'M')
+ ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
+ else // (jobd == 'N')
+ ldwork = 5*l*nobr;
+ }
+ else // meth == 'N' && (batch == 'L' || batch == 'O')
+ {
+ ldwork = 5*(m+l)*nobr + 1;
+ }
+ }
+ else if (alg == 'F')
+ {
+/*
+For the second LDWORK case, code and documentation don't match:
+doc line 276: BATCH = 'F', 'I'
+code line 586: BATCH = 'F', 'I', 'O'
+The third case with BATCH = 'O' is never reached.
+
+
+IB01AD.f Lines 273-279:
+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';
+
+
+IB01AD.f Lines 499-500:
+ ONEBCH = LSAME( BATCH, 'O' )
+ FIRST = LSAME( BATCH, 'F' ) .OR. ONEBCH
+
+
+IB01AD.f Lines 583-591:
+ ELSE IF ( FQRALG ) THEN
+ IF ( .NOT.ONEBCH .AND. CONNEC ) THEN
+ MINWRK = NR*( M + L + 3 )
+ ELSE IF ( FIRST .OR. INTERM ) THEN // (batch = F || O) || batch = I
+ MINWRK = NR*( M + L + 1 ) ^
+ ELSE |
+ MINWRK = 2*NR*( M + L + 1 ) + NR ???
+ END IF
+ ELSE
+*/
+ if (batch != 'O' && conct == 'C')
+ ldwork = (m+l)*2*nobr*(m+l+3);
+ else if (batch == 'F' || batch == 'O' || batch == 'I') // && conct == 'N'
+ ldwork = (m+l)*2*nobr*(m+l+1);
+ else // (batch == 'L' && conct == 'N')
+ ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr;
+ }
+ else // (alg == 'Q')
+ {
+ int ns = nsmp - 2*nobr + 1;
+
+ if (ldr >= ns && batch == 'F')
+ {
+ ldwork = 4*(m+l)*nobr;
+ }
+ else if (ldr >= ns && batch == 'O')
+ {
+ if (meth == 'M')
+ ldwork = max (4*(m+l)*nobr, 5*l*nobr);
+ else // (meth == 'N')
+ ldwork = 5*(m+l)*nobr + 1;
+ }
+ else if (conct == 'C' && (batch == 'I' || batch == 'L'))
+ {
+ ldwork = 4*(nobr+1)*(m+l)*nobr;
+ }
+ else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N')
+ {
+ ldwork = 6*(m+l)*nobr;
+ }
+ }
+
+/*
+IB01AD.f Lines 291-195:
+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.
+
+somehow ldrwrk and ldwork must have been mixed up here
+
+*/
+
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicators
+ int iwarn = 0;
+ int info = 0;
+
+
+ // SLICOT routine IB01AD
+ F77_XFCN (ib01ad, IB01AD,
+ (meth, alg, jobd,
+ batch, conct, ctrl,
+ nobr, m, l,
+ nsmp,
+ u.fortran_vec (), ldu,
+ y.fortran_vec (), ldy,
+ n,
+ r.fortran_vec (), ldr,
+ sv.fortran_vec (),
+ rcond, tol,
+ iwork,
+ dwork, ldwork,
+ iwarn, info));
+
+
+ if (f77_exception_encountered)
+ error ("ident: exception in SLICOT subroutine IB01AD");
+
+ static const char* err_msg[] = {
+ "0: OK",
+ "1: a fast algorithm was requested (ALG = 'C', or 'F') "
+ "in sequential data processing, but it failed; the "
+ "routine can be repeatedly called again using the "
+ "standard QR algorithm",
+ "2: the singular value decomposition (SVD) algorithm did "
+ "not converge"};
+
+ static const char* warn_msg[] = {
+ "0: OK",
+ "1: the number of 100 cycles in sequential data "
+ "processing has been exhausted without signaling "
+ "that the last block of data was get; the cycle "
+ "counter was reinitialized",
+ "2: a fast algorithm was requested (ALG = 'C' or 'F'), "
+ "but it failed, and the QR algorithm was then used "
+ "(non-sequential data processing)",
+ "3: all singular values were exactly zero, hence N = 0 "
+ "(both input and output were identically zero)",
+ "4: the least squares problems with coefficient matrix "
+ "U_f, used for computing the weighted oblique "
+ "projection (for METH = 'N'), have a rank-deficient "
+ "coefficient matrix",
+ "5: the least squares problem with coefficient matrix "
+ "r_1 [6], used for computing the weighted oblique "
+ "projection (for METH = 'N'), has a rank-deficient "
+ "coefficient matrix"};
+
+
+ error_msg ("ident", info, 2, err_msg);
+ warning_msg ("ident", iwarn, 5, warn_msg);
+
+
+ // resize
+ int rs = 2*(m+l)*nobr;
+ r.resize (rs, rs);
+
+//////////////////////////////////////////////////////////
+// SLICOT IB01BD - A, B, C, D //
+//////////////////////////////////////////////////////////
+
+ // arguments in
+ char job = 'A';
+ char jobck = 'K';
+
+ // TODO: if meth == 'C', which meth should be taken for IB01AD.f, 'M' or 'N'?
+
+ int nsmpl = nsmp;
+
+ // arguments out
+
+ // error indicators
+ int iwarn = 0;
+ int info = 0;
+
+
+ // SLICOT routine IB01BD
+ F77_XFCN (ib01bd, IB01BD,
+ (meth, job, jobck,
+ nobr, n, m, l,
+ nsmpl,
+ r.fortran_vec (), ldr,
+ a.fortran_vec (), lda,
+ c.fortran_vec (), ldc,
+ b.fortran_vec (), ldb,
+ d.fortran_vec (), ldd,
+ q.fortran_vec (), ldq,
+ ry.fortran_vec (), ldry,
+ s.fortran_vec (), lds,
+ k.fortran_vec (), ldk,
+ tol,
+ iwork,
+ dwork, ldwork,
+ bwork,
+ iwarn, info));
+
+
+ if (f77_exception_encountered)
+ error ("ident: exception in SLICOT subroutine IB01BD");
+
+ static const char* err_msg_b[] = {
+ "0: OK",
+ "1: error message not specified",
+ "2: the singular value decomposition (SVD) algorithm did "
+ "not converge",
+ "3: a singular upper triangular matrix was found",
+ "4: matrix A is (numerically) singular in discrete-"
+ "time case",
+ "5: the Hamiltonian or symplectic matrix H cannot be "
+ "reduced to real Schur form",
+ "6: the real Schur form of the Hamiltonian or "
+ "symplectic matrix H cannot be appropriately ordered",
+ "7: the Hamiltonian or symplectic matrix H has less "
+ "than N stable eigenvalues",
+ "8: the N-th order system of linear algebraic "
+ "equations, from which the solution matrix X would "
+ "be obtained, is singular to working precision",
+ "9: the QR algorithm failed to complete the reduction "
+ "of the matrix Ac to Schur canonical form, T",
+ "10: the QR algorithm did not converge"};
+
+ static const char* warn_msg_b[] = {
+ "0: OK",
+ "1: warning message not specified",
+ "2: warning message not specified",
+ "3: warning message not specified",
+ "4: a least squares problem to be solved has a "
+ "rank-deficient coefficient matrix",
+ "5: the computed covariance matrices are too small. "
+ "The problem seems to be a deterministic one; the "
+ "gain matrix is set to zero"};
+
+
+ error_msg ("ident", info, 10, err_msg_b);
+ warning_msg ("ident", iwarn, 5, warn_msg_b);
+
+
+ // return values
+ //retval(0) = octave_value (n);
+ //retval(1) = r;
+ //retval(2) = sv;
+ }
+
+ return retval;
+}