Menu
▾
▴
octave-cvsupdate
|
From: <par...@us...> - 2012-03-15 20:04:56
|
Revision: 9906
http://octave.svn.sourceforge.net/octave/?rev=9906&view=rev
Author: paramaniac
Date: 2012-03-15 20:04:50 +0000 (Thu, 15 Mar 2012)
Log Message:
-----------
control-devel: quicksave id draft code (2)
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/src/slident.cc
Modified: trunk/octave-forge/extra/control-devel/src/slident.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-15 19:24:20 UTC (rev 9905)
+++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-15 20:04:50 UTC (rev 9906)
@@ -79,7 +79,7 @@
int nargin = args.length ();
octave_value_list retval;
- if (nargin != 11)
+ if (nargin != 12)
{
print_usage ();
}
@@ -106,6 +106,7 @@
double rcond = args(9).double_value ();
double tol = args(10).double_value ();
+ double tolb = args(11).double_value ();
if (imeth == 0)
@@ -389,7 +390,35 @@
int nsmpl = nsmp;
// arguments out
+ lda = max (1, n);
+ ldc = max (1, l);
+ ldb = max (1, n);
+ ldd = max (1, l);
+ ldq = n; // if JOBCK = 'C' or 'K'
+ ldry = l; // if JOBCK = 'C' or 'K'
+ lds = n; // if JOBCK = 'C' or 'K'
+ ldk = n; // if JOBCK = 'K'
+
+ Matrix a (lda, n);
+ Matrix c (ldc, n);
+ Matrix b (ldb, m);
+ Matrix d (ldd, m);
+
+ Matrix q (ldq, n);
+ Matrix ry (ldry, l);
+ Matrix s (lds, l);
+ Matrix k (ldk, l);
+
+ // workspace
+ int liwork;
+ int ldwork;
+
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+
// error indicators
int iwarn = 0;
int info = 0;
@@ -409,7 +438,7 @@
ry.fortran_vec (), ldry,
s.fortran_vec (), lds,
k.fortran_vec (), ldk,
- tol,
+ tolb,
iwork,
dwork, ldwork,
bwork,
@@ -455,8 +484,27 @@
error_msg ("ident", info, 10, err_msg_b);
warning_msg ("ident", iwarn, 5, warn_msg_b);
+ // resize
+ a.resize (n, n);
+ c.resize (l, n);
+ b.resize (n, m);
+ d.resize (l, m);
+ q.resize (n, n);
+ ry.resize (l, l);
+ s.resize (n, l);
+ k.resize (n, l);
+
// return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = c;
+ retval(3) = d;
+
+ retval(4) = q;
+ retval(5) = ry;
+ retval(6) = s;
+ retval(7) = k;
//retval(0) = octave_value (n);
//retval(1) = r;
//retval(2) = sv;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-03-17 07:20:57
|
Revision: 9930
http://octave.svn.sourceforge.net/octave/?rev=9930&view=rev
Author: paramaniac
Date: 2012-03-17 07:20:50 +0000 (Sat, 17 Mar 2012)
Log Message:
-----------
control-devel: prepare computation of initial state
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/src/slident.cc
Modified: trunk/octave-forge/extra/control-devel/src/slident.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-17 02:55:09 UTC (rev 9929)
+++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-17 07:20:50 UTC (rev 9930)
@@ -67,6 +67,22 @@
bool* BWORK,
int& IWARN, int& INFO);
+ int F77_FUNC (ib01cd, IB01CD)
+ (char& JOBX0, char& COMUSE, char& JOB,
+ int& N, int& M, int& L,
+ int& NSMP,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* D, int& LDD,
+ double* U, int& LDU,
+ double* Y, int& LDY,
+ double* X0,
+ double* V, int& LDV,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
}
// PKG_ADD: autoload ("slident", "devel_slicot_functions.oct");
@@ -85,6 +101,10 @@
}
else
{
+////////////////////////////////////////////////////////////////////////////////////
+// SLICOT IB01AD - preprocess the input-output data //
+////////////////////////////////////////////////////////////////////////////////////
+
// arguments in
char meth;
char alg;
@@ -359,9 +379,9 @@
int rs = 2*(m+l)*nobr;
r.resize (rs, rs);
-//////////////////////////////////////////////////////////
-// SLICOT IB01BD - A, B, C, D //
-//////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////
+// SLICOT IB01BD - estimating system matrices, Kalman gain, and covariances //
+////////////////////////////////////////////////////////////////////////////////////
// arguments in
char job = 'A';
@@ -532,6 +552,57 @@
ry.resize (l, l);
s.resize (n, l);
k.resize (n, l);
+
+////////////////////////////////////////////////////////////////////////////////////
+// SLICOT IB01CD - estimating the initial state //
+////////////////////////////////////////////////////////////////////////////////////
+
+ // SLICOT routine IB01CD
+ F77_XFCN (ib01cd, IB01CD,
+ (jobx0, comuse, job***,
+ n, m, l,
+ nsmp,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ c.fortran_vec (), ldc,
+ d.fortran_vec (), ldd,
+ u.fortran_vec (), ldu,
+ y.fortran_vec (), ldy,
+ x0.fortran_vec (),
+ v.fortran_vec (), ldv,
+ tolb,
+ iwork_c,
+ dwork_c, ldwork_c,
+ iwarn_c, info_c));
+
+
+ if (f77_exception_encountered)
+ error ("ident: exception in SLICOT subroutine IB01CD");
+
+ static const char* err_msg_c[] = {
+ "0: OK",
+ "1: the QR algorithm failed to compute all the "
+ "eigenvalues of the matrix A (see LAPACK Library "
+ "routine DGEES); the locations DWORK(i), for "
+ "i = g+1:g+N*N, contain the partially converged "
+ "Schur form",
+ "2: the singular value decomposition (SVD) algorithm did "
+ "not converge"};
+
+ static const char* warn_msg_c[] = {
+ "0: OK",
+ "1: warning message not specified",
+ "2: warning message not specified",
+ "3: warning message not specified",
+ "4: the least squares problem to be solved has a "
+ "rank-deficient coefficient matrix",
+ "5: warning message not specified",
+ "6: the matrix A is unstable; the estimated x(0) "
+ "and/or B and D could be inaccurate"};
+
+
+ error_msg ("ident", info_c, 2, err_msg_c);
+ warning_msg ("ident", iwarn_c, 6, warn_msg_c);
// return values
retval(0) = a;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
[Octave-cvsupdate] SF.net SVN: octave:[10423]
trunk/octave-forge/extra/control-devel/src/ slident.cc
From: <par...@us...> - 2012-05-12 14:02:27
|
Revision: 10423
http://octave.svn.sourceforge.net/octave/?rev=10423&view=rev
Author: paramaniac
Date: 2012-05-12 14:02:19 +0000 (Sat, 12 May 2012)
Log Message:
-----------
control-devel: clean up sliding further
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/src/slident.cc
Modified: trunk/octave-forge/extra/control-devel/src/slident.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-12 13:50:09 UTC (rev 10422)
+++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-12 14:02:19 UTC (rev 10423)
@@ -329,7 +329,7 @@
// warning ("==================== ldwork_a before: %d =====================", ldwork_a);
// ldwork_a = (ns+2)*(2*(m+l)*nobr);
-ldwork_a = max (ldwork_a, (ns+2)*(2*(m+l)*nobr));
+//////////ldwork_a = max (ldwork_a, (ns+2)*(2*(m+l)*nobr));
// ldwork_a *= 3;
// warning ("==================== ldwork_a after: %d =====================", ldwork_a);
@@ -470,7 +470,7 @@
int ldw1;
int ldw2;
int ldw3;
-/*
+
if (meth_b == 'M')
{
int ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
@@ -514,57 +514,10 @@
ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1);
}
-*/
-
- int ldw1ax = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
- int ldw1bx = max (2*(l*nobr-l)*n+n*n+7*n,
- (l*nobr-l)*n+n+6*m*nobr,
- (l*nobr-l)*n+n+max (l+m*nobr, l*nobr + max (3*l*nobr+1, m)));
- int ldw1x = max (ldw1ax, ldw1bx);
- int aw;
-
- if (m == 0 || job == 'C')
- aw = n + n*n;
- else
- aw = 0;
-
- int ldw2x = l*nobr*n + max ((l*nobr-l)*n+aw+2*n+max(5*n,(2*m+l)*nobr+l), 4*(m*nobr+n)+1, m*nobr+2*n+l );
-
-
-
- int ldw1y = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l,
- 2*(l*nobr-l)*n+n*n+8*n,
- n+4*(m*nobr+n)+1,
- m*nobr+3*n+l);
- int ldw2y;
- if (m == 0 || job == 'C')
- int ldw2y = 0;
- else
- int ldw2y = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1);
-
-
- int ldw1az = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
- int ldw1bz = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l,
- 2*(l*nobr-l)*n+n*n+8*n,
- n+4*(m*nobr+n)+1,
- m*nobr+3*n+l);
-
- int ldw1z = max (ldw1az, ldw1bz);
-
- int ldw2z = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1);
-
-
- ldw1 = max (ldw1x, ldw1y, ldw1z);
- ldw2 = max (ldw2x, ldw2y, ldw2z);
-
-
-
ldw3 = max(4*n*n + 2*n*l + l*l + max (3*l, n*l), 14*n*n + 12*n + 5);
ldwork_b = max (ldw1, ldw2, ldw3);
- //
- ldwork_b *= 3;
OCTAVE_LOCAL_BUFFER (int, iwork_b, liwork_b);
OCTAVE_LOCAL_BUFFER (double, dwork_b, ldwork_b);
@@ -745,10 +698,6 @@
retval(7) = k;
retval(8) = x0;
- //retval(8) = ColumnVector (n);
- //retval(0) = octave_value (n);
- //retval(1) = r;
- //retval(2) = sv;
}
return retval;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
[Octave-cvsupdate] SF.net SVN: octave:[10457]
trunk/octave-forge/extra/control-devel/src/ slident.cc
From: <par...@us...> - 2012-05-18 16:04:27
|
Revision: 10457
http://octave.svn.sourceforge.net/octave/?rev=10457&view=rev
Author: paramaniac
Date: 2012-05-18 16:04:20 +0000 (Fri, 18 May 2012)
Log Message:
-----------
control-devel: work on multi-experiment identification
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/src/slident.cc
Modified: trunk/octave-forge/extra/control-devel/src/slident.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-18 15:01:56 UTC (rev 10456)
+++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-18 16:04:20 UTC (rev 10457)
@@ -28,7 +28,8 @@
*/
#include <octave/oct.h>
-#include <f77-fcn.h>
+#include <octave/f77-fcn.h>
+#include <octave/Cell.h>
#include "common.h"
extern "C"
@@ -114,14 +115,16 @@
char conct;
char ctrl;
- Matrix y = args(0).matrix_value ();
- Matrix u = args(1).matrix_value ();
+ const Cell y_cell = args(0).cell_value ();
+ const Cell u_cell = args(1).cell_value ();
+ //Matrix y = args(0).matrix_value ();
+ //Matrix u = args(1).matrix_value ();
int nobr = args(2).int_value ();
int nuser = args(3).int_value ();
const int imeth = args(4).int_value ();
const int ialg = args(5).int_value ();
- const int ibatch = args(6).int_value ();
+ // const int ibatch = args(6).int_value (); // löschen
const int iconct = args(7).int_value ();
const int ictrl = args(8).int_value ();
@@ -169,24 +172,6 @@
jobd = 'M';
else // meth_a == 'N'
jobd = 'N'; // IB01AD.f says: This parameter is not relevant for METH = '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';
@@ -199,32 +184,66 @@
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 n_exp = y_cell.nelem (); // number of experiments
- if (batch == 'O')
- {
- if (nsmp < 2*(m+l+1)*nobr - 1)
- error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1");
- }
+ int m = u_cell.elem(0).columns (); // m: number of inputs
+ int l = y_cell.elem(0).columns (); // l: number of outputs
+
+ // arguments out
+ int n;
+ int ldr;
+
+ if (meth_a == 'M' && jobd == 'M')
+ ldr = max (2*(m+l)*nobr, 3*m*nobr);
+ else if (meth_a == 'N' || (meth_a == '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);
+
+
+ for (int i = 0; i < n_exp; i++)
{
- if (nsmp < 2*nobr)
- error ("slident: require NSMP >= 2*NOBR");
- }
+ if (n_exp == 1)
+ batch = 'O'; // one block only
+ else if (i == 0)
+ batch = 'F'; // first block
+ else if (i == n_exp-1)
+ batch = 'L'; // last block
+ else
+ batch = 'I'; // intermediate block
+
+ Matrix y = y_cell.elem(i).matrix_value ();
+ Matrix u = u_cell.elem(i).matrix_value ();
+
+ //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 (batch == 'O')
+ {
+ if (nsmp < 2*(m+l+1)*nobr - 1)
+ error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1");
+ }
+ else
+ {
+ if (nsmp < 2*nobr)
+ error ("slident: require NSMP >= 2*NOBR");
+ }
- if (m == 0)
- ldu = 1;
- else // m > 0
- ldu = nsmp;
+ int ldu;
+
+ if (m == 0)
+ ldu = 1;
+ else // m > 0
+ ldu = nsmp;
- int ldy = nsmp;
-
+ int ldy = nsmp;
+/*
// arguments out
int n;
int ldr;
@@ -238,78 +257,78 @@
Matrix r (ldr, 2*(m+l)*nobr);
ColumnVector sv (l*nobr);
+*/
+ // workspace
+ int liwork_a;
- // workspace
- int liwork_a;
+ if (meth_a == 'N') // if METH = 'N'
+ liwork_a = (m+l)*nobr;
+ else if (alg == 'F') // if METH = 'M' and ALG = 'F'
+ liwork_a = m+l;
+ else // if METH = 'M' and ALG = 'C' or 'Q'
+ liwork_a = 0;
- if (meth_a == 'N') // if METH = 'N'
- liwork_a = (m+l)*nobr;
- else if (alg == 'F') // if METH = 'M' and ALG = 'F'
- liwork_a = m+l;
- else // if METH = 'M' and ALG = 'C' or 'Q'
- liwork_a = 0;
+ // TODO: Handle 'k' for DWORK
- // TODO: Handle 'k' for DWORK
-
- int ldwork_a;
- int ns = nsmp - 2*nobr + 1;
+ int ldwork_a;
+ //int ns = nsmp - 2*nobr + 1;
- if (alg == 'C')
- {
- if (batch == 'F' || batch == 'I')
+ if (alg == 'C')
{
- if (conct == 'C')
- ldwork_a = (4*nobr-2)*(m+l);
- else // (conct == 'N')
- ldwork_a = 1;
- }
- else if (meth_a == 'M') // && (batch == 'L' || batch == 'O')
- {
- if (conct == 'C' && batch == 'L')
- ldwork_a = max ((4*nobr-2)*(m+l), 5*l*nobr);
- else if (jobd == 'M')
- ldwork_a = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
- else // (jobd == 'N')
- ldwork_a = 5*l*nobr;
- }
- else // meth_b == 'N' && (batch == 'L' || batch == 'O')
- {
- ldwork_a = 5*(m+l)*nobr + 1;
- }
- }
- else if (alg == 'F')
- {
- if (batch != 'O' && conct == 'C')
- ldwork_a = (m+l)*2*nobr*(m+l+3);
- else if (batch == 'F' || batch == 'I') // && conct == 'N'
- ldwork_a = (m+l)*2*nobr*(m+l+1);
- else // (batch == 'L' || '0' && conct == 'N')
- ldwork_a = (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_a = 4*(m+l)*nobr;
- }
- else if (ldr >= ns && batch == 'O')
- {
- if (meth_a == 'M')
- ldwork_a = max (4*(m+l)*nobr, 5*l*nobr);
- else // (meth == 'N')
+ if (batch == 'F' || batch == 'I')
+ {
+ if (conct == 'C')
+ ldwork_a = (4*nobr-2)*(m+l);
+ else // (conct == 'N')
+ ldwork_a = 1;
+ }
+ else if (meth_a == 'M') // && (batch == 'L' || batch == 'O')
+ {
+ if (conct == 'C' && batch == 'L')
+ ldwork_a = max ((4*nobr-2)*(m+l), 5*l*nobr);
+ else if (jobd == 'M')
+ ldwork_a = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
+ else // (jobd == 'N')
+ ldwork_a = 5*l*nobr;
+ }
+ else // meth_b == 'N' && (batch == 'L' || batch == 'O')
+ {
ldwork_a = 5*(m+l)*nobr + 1;
+ }
}
- else if (conct == 'C' && (batch == 'I' || batch == 'L'))
+ else if (alg == 'F')
{
- ldwork_a = 4*(nobr+1)*(m+l)*nobr;
+ if (batch != 'O' && conct == 'C')
+ ldwork_a = (m+l)*2*nobr*(m+l+3);
+ else if (batch == 'F' || batch == 'I') // && conct == 'N'
+ ldwork_a = (m+l)*2*nobr*(m+l+1);
+ else // (batch == 'L' || '0' && conct == 'N')
+ ldwork_a = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr;
}
- else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N')
+ else // (alg == 'Q')
{
- ldwork_a = 6*(m+l)*nobr;
+ int ns = nsmp - 2*nobr + 1;
+
+ if (ldr >= ns && batch == 'F')
+ {
+ ldwork_a = 4*(m+l)*nobr;
+ }
+ else if (ldr >= ns && batch == 'O')
+ {
+ if (meth_a == 'M')
+ ldwork_a = max (4*(m+l)*nobr, 5*l*nobr);
+ else // (meth == 'N')
+ ldwork_a = 5*(m+l)*nobr + 1;
+ }
+ else if (conct == 'C' && (batch == 'I' || batch == 'L'))
+ {
+ ldwork_a = 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_a = 6*(m+l)*nobr;
+ }
}
- }
/*
IB01AD.f Lines 438-445
@@ -343,66 +362,67 @@
*/
- OCTAVE_LOCAL_BUFFER (int, iwork_a, liwork_a);
- OCTAVE_LOCAL_BUFFER (double, dwork_a, ldwork_a);
+ OCTAVE_LOCAL_BUFFER (int, iwork_a, liwork_a);
+ OCTAVE_LOCAL_BUFFER (double, dwork_a, ldwork_a);
- // error indicators
- int iwarn_a = 0;
- int info_a = 0;
+ // error indicators
+ int iwarn_a = 0;
+ int info_a = 0;
- // SLICOT routine IB01AD
- F77_XFCN (ib01ad, IB01AD,
- (meth_a, 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_a,
- iwork_a,
- dwork_a, ldwork_a,
- iwarn_a, info_a));
+ // SLICOT routine IB01AD
+ F77_XFCN (ib01ad, IB01AD,
+ (meth_a, 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_a,
+ iwork_a,
+ dwork_a, ldwork_a,
+ iwarn_a, info_a));
- if (f77_exception_encountered)
- error ("ident: exception in SLICOT subroutine IB01AD");
+ 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* 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"};
+ 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: IB01AD", info_a, 2, err_msg);
- warning_msg ("ident: IB01AD", iwarn_a, 5, warn_msg);
+ error_msg ("ident: IB01AD", info_a, 2, err_msg);
+ warning_msg ("ident: IB01AD", iwarn_a, 5, warn_msg);
+ }
// resize
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
[Octave-cvsupdate] SF.net SVN: octave:[10418]
trunk/octave-forge/extra/control-devel/src/ slident.cc
From: <par...@us...> - 2012-05-12 13:07:23
|
Revision: 10418
http://octave.svn.sourceforge.net/octave/?rev=10418&view=rev
Author: paramaniac
Date: 2012-05-12 13:07:17 +0000 (Sat, 12 May 2012)
Log Message:
-----------
control-devel: clean up variable names
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/src/slident.cc
Modified: trunk/octave-forge/extra/control-devel/src/slident.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-12 13:02:16 UTC (rev 10417)
+++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-12 13:07:17 UTC (rev 10418)
@@ -106,14 +106,13 @@
////////////////////////////////////////////////////////////////////////////////////
// arguments in
- char meth;
+ char meth_a;
+ char meth_b;
char alg;
char jobd;
char batch;
char conct;
char ctrl;
- char metha;
- char jobda; // ??? unused
Matrix y = args(0).matrix_value ();
Matrix u = args(1).matrix_value ();
@@ -128,23 +127,26 @@
const int ictrl = args(9).int_value ();
double rcond = args(10).double_value ();
- double tol = args(11).double_value ();
- double tolb = args(10).double_value (); // tolb = rcond
+ double tol_a = args(11).double_value ();
+ double tol_b = rcond;
+ double tol_c = rcond;
+ //double tol_b = args(10).double_value (); // tol_b = rcond
+
switch (imeth)
{
case 0:
- meth = 'M';
- metha = 'M';
+ meth_b = 'M';
+ meth_a = 'M';
break;
case 1:
- meth = 'N';
- metha = 'N';
+ meth_b = 'N';
+ meth_a = 'N';
break;
case 2:
- meth = 'C';
- metha = 'N'; // no typo here
+ meth_b = 'C';
+ meth_a = 'N'; // no typo here
break;
default:
error ("slib01ad: argument 'meth' invalid");
@@ -165,7 +167,7 @@
error ("slib01ad: argument 'alg' invalid");
}
- if (meth == 'C')
+ if (meth_b == 'C')
jobd = 'N';
else if (ijobd == 0)
jobd = 'M';
@@ -231,9 +233,9 @@
int n;
int ldr;
- if (metha == 'M' && jobd == 'M')
+ if (meth_a == 'M' && jobd == 'M')
ldr = max (2*(m+l)*nobr, 3*m*nobr);
- else if (metha == 'N' || (metha == 'M' && jobd == 'N'))
+ else if (meth_a == 'N' || (meth_a == 'M' && jobd == 'N'))
ldr = 2*(m+l)*nobr;
else
error ("slib01ad: could not handle 'ldr' case");
@@ -242,18 +244,18 @@
ColumnVector sv (l*nobr);
// workspace
- int liwork;
+ int liwork_a;
- if (metha == 'N') // if METH = 'N'
- liwork = (m+l)*nobr;
+ if (meth_a == 'N') // if METH = 'N'
+ liwork_a = (m+l)*nobr;
else if (alg == 'F') // if METH = 'M' and ALG = 'F'
- liwork = m+l;
+ liwork_a = m+l;
else // if METH = 'M' and ALG = 'C' or 'Q'
- liwork = 0;
+ liwork_a = 0;
// TODO: Handle 'k' for DWORK
- int ldwork;
+ int ldwork_a;
int ns = nsmp - 2*nobr + 1;
if (alg == 'C')
@@ -261,32 +263,32 @@
if (batch == 'F' || batch == 'I')
{
if (conct == 'C')
- ldwork = (4*nobr-2)*(m+l);
+ ldwork_a = (4*nobr-2)*(m+l);
else // (conct == 'N')
- ldwork = 1;
+ ldwork_a = 1;
}
- else if (metha == 'M') // && (batch == 'L' || batch == 'O')
+ else if (meth_a == 'M') // && (batch == 'L' || batch == 'O')
{
if (conct == 'C' && batch == 'L')
- ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr);
+ ldwork_a = 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);
+ ldwork_a = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
else // (jobd == 'N')
- ldwork = 5*l*nobr;
+ ldwork_a = 5*l*nobr;
}
- else // meth == 'N' && (batch == 'L' || batch == 'O')
+ else // meth_b == 'N' && (batch == 'L' || batch == 'O')
{
- ldwork = 5*(m+l)*nobr + 1;
+ ldwork_a = 5*(m+l)*nobr + 1;
}
}
else if (alg == 'F')
{
if (batch != 'O' && conct == 'C')
- ldwork = (m+l)*2*nobr*(m+l+3);
+ ldwork_a = (m+l)*2*nobr*(m+l+3);
else if (batch == 'F' || batch == 'I') // && conct == 'N'
- ldwork = (m+l)*2*nobr*(m+l+1);
+ ldwork_a = (m+l)*2*nobr*(m+l+1);
else // (batch == 'L' || '0' && conct == 'N')
- ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr;
+ ldwork_a = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr;
}
else // (alg == 'Q')
{
@@ -294,22 +296,22 @@
if (ldr >= ns && batch == 'F')
{
- ldwork = 4*(m+l)*nobr;
+ ldwork_a = 4*(m+l)*nobr;
}
else if (ldr >= ns && batch == 'O')
{
- if (metha == 'M')
- ldwork = max (4*(m+l)*nobr, 5*l*nobr);
+ if (meth_a == 'M')
+ ldwork_a = max (4*(m+l)*nobr, 5*l*nobr);
else // (meth == 'N')
- ldwork = 5*(m+l)*nobr + 1;
+ ldwork_a = 5*(m+l)*nobr + 1;
}
else if (conct == 'C' && (batch == 'I' || batch == 'L'))
{
- ldwork = 4*(nobr+1)*(m+l)*nobr;
+ ldwork_a = 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;
+ ldwork_a = 6*(m+l)*nobr;
}
}
@@ -325,11 +327,11 @@
C cache size is large enough to accommodate R, U, Y, and DWORK.
*/
-// warning ("==================== ldwork before: %d =====================", ldwork);
-// ldwork = (ns+2)*(2*(m+l)*nobr);
-ldwork = max (ldwork, (ns+2)*(2*(m+l)*nobr));
-// ldwork *= 3;
-// warning ("==================== ldwork after: %d =====================", ldwork);
+// warning ("==================== ldwork_a before: %d =====================", ldwork_a);
+// ldwork_a = (ns+2)*(2*(m+l)*nobr);
+ldwork_a = max (ldwork_a, (ns+2)*(2*(m+l)*nobr));
+// ldwork_a *= 3;
+// warning ("==================== ldwork_a after: %d =====================", ldwork_a);
/*
@@ -345,17 +347,17 @@
*/
- OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
- OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+ OCTAVE_LOCAL_BUFFER (int, iwork_a, liwork_a);
+ OCTAVE_LOCAL_BUFFER (double, dwork_a, ldwork_a);
// error indicators
- int iwarn = 0;
- int info = 0;
+ int iwarn_a = 0;
+ int info_a = 0;
// SLICOT routine IB01AD
F77_XFCN (ib01ad, IB01AD,
- (metha, alg, jobd,
+ (meth_a, alg, jobd,
batch, conct, ctrl,
nobr, m, l,
nsmp,
@@ -364,10 +366,10 @@
n,
r.fortran_vec (), ldr,
sv.fortran_vec (),
- rcond, tol,
- iwork,
- dwork, ldwork,
- iwarn, info));
+ rcond, tol_a,
+ iwork_a,
+ dwork_a, ldwork_a,
+ iwarn_a, info_a));
if (f77_exception_encountered)
@@ -403,8 +405,8 @@
"coefficient matrix"};
- error_msg ("ident", info, 2, err_msg);
- warning_msg ("ident", iwarn, 5, warn_msg);
+ error_msg ("ident", info_a, 2, err_msg);
+ warning_msg ("ident", iwarn_a, 5, warn_msg);
// resize
@@ -430,8 +432,6 @@
char job = 'A';
char jobck = 'K';
- // TODO: if meth == 'C', which meth should be taken for IB01AD.f, 'M' or 'N'?
-
int nsmpl = nsmp;
if (nsmpl < 2*(m+l)*nobr)
@@ -471,7 +471,7 @@
int ldw2;
int ldw3;
/*
- if (meth == 'M')
+ if (meth_b == 'M')
{
int ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
int ldw1b = max (2*(l*nobr-l)*n+n*n+7*n,
@@ -488,7 +488,7 @@
ldw2 = l*nobr*n + max ((l*nobr-l)*n+aw+2*n+max(5*n,(2*m+l)*nobr+l), 4*(m*nobr+n)+1, m*nobr+2*n+l );
}
- else if (meth == 'N')
+ else if (meth_b == 'N')
{
ldw1 = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l,
2*(l*nobr-l)*n+n*n+8*n,
@@ -501,7 +501,7 @@
ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1);
}
- else // (meth == 'C')
+ else // (meth_b == 'C')
{
int ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
int ldw1b = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l,
@@ -578,7 +578,7 @@
// SLICOT routine IB01BD
F77_XFCN (ib01bd, IB01BD,
- (meth, job, jobck,
+ (meth_b, job, jobck,
nobr, n, m, l,
nsmpl,
r.fortran_vec (), ldr,
@@ -590,7 +590,7 @@
ry.fortran_vec (), ldry,
s.fortran_vec (), lds,
k.fortran_vec (), ldk,
- tolb,
+ tol_b,
iwork_b,
dwork_b, ldwork_b,
bwork,
@@ -698,7 +698,7 @@
y.fortran_vec (), ldy,
x0.fortran_vec (),
v.fortran_vec (), ldv,
- tolb,
+ tol_c,
iwork_c,
dwork_c, ldwork_c,
iwarn_c, info_c));
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
[Octave-cvsupdate] SF.net SVN: octave:[10488]
trunk/octave-forge/extra/control-devel/src/ slident.cc
From: <par...@us...> - 2012-05-22 06:31:56
|
Revision: 10488
http://octave.svn.sourceforge.net/octave/?rev=10488&view=rev
Author: paramaniac
Date: 2012-05-22 06:31:46 +0000 (Tue, 22 May 2012)
Log Message:
-----------
control-devel: increase workspace for better efficiency, add comments
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/src/slident.cc
Modified: trunk/octave-forge/extra/control-devel/src/slident.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-22 05:59:21 UTC (rev 10487)
+++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-22 06:31:46 UTC (rev 10488)
@@ -180,7 +180,7 @@
else
ctrl = 'N';
-
+ // m and l are equal for all experiments, checked by iddata class
int n_exp = y_cell.nelem (); // number of experiments
int m = u_cell.elem(0).columns (); // m: number of inputs
int l = y_cell.elem(0).columns (); // l: number of outputs
@@ -201,6 +201,7 @@
ColumnVector sv (l*nobr);
+ // repeat for every experiment in the dataset
for (int i = 0; i < n_exp; i++)
{
if (n_exp == 1)
@@ -215,13 +216,13 @@
Matrix y = y_cell.elem(i).matrix_value ();
Matrix u = u_cell.elem(i).matrix_value ();
- //int m = u.columns (); // m: number of inputs
- //int l = y.columns (); // l: number of outputs
- int nsmp = y.rows (); // nsmp: number of samples
- nsmpl += nsmp;
// y.rows == u.rows is checked by iddata class
- // TODO: check minimal nsmp size
-
+ // int m = u.columns (); // m: number of inputs
+ // int l = y.columns (); // l: number of outputs
+ int nsmp = y.rows (); // nsmp: number of samples in the current experiment
+ nsmpl += nsmp; // nsmpl: total number of samples of all experiments
+
+ // minimal nsmp size checked by __slicot_identification__.m
if (batch == 'O')
{
if (nsmp < 2*(m+l+1)*nobr - 1)
@@ -255,6 +256,7 @@
// TODO: Handle 'k' for DWORK
int ldwork_a;
+ int ns = nsmp - 2*nobr + 1;
if (alg == 'C')
{
@@ -290,7 +292,7 @@
}
else // (alg == 'Q')
{
- int ns = nsmp - 2*nobr + 1;
+ // int ns = nsmp - 2*nobr + 1;
if (ldr >= ns && batch == 'F')
{
@@ -313,38 +315,33 @@
}
}
-/*
-IB01AD.f Lines 438-445
-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.
-*/
+ /*
+ IB01AD.f Lines 438-445
+ 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.
+ */
-// warning ("==================== ldwork_a before: %d =====================", ldwork_a);
-// ldwork_a = (ns+2)*(2*(m+l)*nobr);
-//////////ldwork_a = max (ldwork_a, (ns+2)*(2*(m+l)*nobr));
-// ldwork_a *= 3;
-// warning ("==================== ldwork_a after: %d =====================", ldwork_a);
+ ldwork_a = max (ldwork_a, (ns+2)*(2*(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.
-/*
-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
-somehow ldrwrk and ldwork must have been mixed up here
+ */
-*/
-
OCTAVE_LOCAL_BUFFER (int, iwork_a, liwork_a);
OCTAVE_LOCAL_BUFFER (double, dwork_a, ldwork_a);
@@ -609,8 +606,10 @@
char jobbd = 'D';
// arguments out
- Cell x0_cell (n_exp, 1);
-
+ Cell x0_cell (n_exp, 1); // cell of initial state vectors x0
+
+ // repeat for every experiment in the dataset
+ // compute individual initial state vector x0 for every experiment
for (int i = 0; i < n_exp; i++)
{
Matrix y = y_cell.elem(i).matrix_value ();
@@ -628,7 +627,7 @@
int ldy = nsmp;
-
+ // arguments out
ColumnVector x0 (n);
Matrix v (ldv, n);
@@ -697,7 +696,7 @@
error_msg ("ident: IB01CD", info_c, 2, err_msg_c);
warning_msg ("ident: IB01CD", iwarn_c, 6, warn_msg_c);
- x0_cell.elem(i) = x0;
+ x0_cell.elem(i) = x0; // add x0 from the current experiment to cell of initial state vectors
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X