Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10988] trunk/octave-forge/main/control/src
|
From: <par...@us...> - 2012-09-09 19:09:57
|
Revision: 10988
http://octave.svn.sourceforge.net/octave/?rev=10988&view=rev
Author: paramaniac
Date: 2012-09-09 19:09:48 +0000 (Sun, 09 Sep 2012)
Log Message:
-----------
control: add underscores to C++ filenames
Modified Paths:
--------------
trunk/octave-forge/main/control/src/__control_slicot_functions__.cc
Added Paths:
-----------
trunk/octave-forge/main/control/src/sl_ab01od.cc
trunk/octave-forge/main/control/src/sl_ab04md.cc
trunk/octave-forge/main/control/src/sl_ab08nd.cc
trunk/octave-forge/main/control/src/sl_ab09hd.cc
trunk/octave-forge/main/control/src/sl_ab09id.cc
trunk/octave-forge/main/control/src/sl_ab09jd.cc
trunk/octave-forge/main/control/src/sl_ab13ad.cc
trunk/octave-forge/main/control/src/sl_ab13bd.cc
trunk/octave-forge/main/control/src/sl_ab13dd.cc
trunk/octave-forge/main/control/src/sl_ag08bd.cc
trunk/octave-forge/main/control/src/sl_ib01ad.cc
trunk/octave-forge/main/control/src/sl_ib01cd.cc
trunk/octave-forge/main/control/src/sl_ident.cc
trunk/octave-forge/main/control/src/sl_sb01bd.cc
trunk/octave-forge/main/control/src/sl_sb02od.cc
trunk/octave-forge/main/control/src/sl_sb03md.cc
trunk/octave-forge/main/control/src/sl_sb03od.cc
trunk/octave-forge/main/control/src/sl_sb04md.cc
trunk/octave-forge/main/control/src/sl_sb04qd.cc
trunk/octave-forge/main/control/src/sl_sb10dd.cc
trunk/octave-forge/main/control/src/sl_sb10ed.cc
trunk/octave-forge/main/control/src/sl_sb10fd.cc
trunk/octave-forge/main/control/src/sl_sb10hd.cc
trunk/octave-forge/main/control/src/sl_sb10id.cc
trunk/octave-forge/main/control/src/sl_sb10jd.cc
trunk/octave-forge/main/control/src/sl_sb10kd.cc
trunk/octave-forge/main/control/src/sl_sb10yd.cc
trunk/octave-forge/main/control/src/sl_sb10zd.cc
trunk/octave-forge/main/control/src/sl_sb16ad.cc
trunk/octave-forge/main/control/src/sl_sb16bd.cc
trunk/octave-forge/main/control/src/sl_sb16cd.cc
trunk/octave-forge/main/control/src/sl_sg02ad.cc
trunk/octave-forge/main/control/src/sl_sg03ad.cc
trunk/octave-forge/main/control/src/sl_sg03bd.cc
trunk/octave-forge/main/control/src/sl_tb01id.cc
trunk/octave-forge/main/control/src/sl_tb01pd.cc
trunk/octave-forge/main/control/src/sl_tb01ud.cc
trunk/octave-forge/main/control/src/sl_tb04bd.cc
trunk/octave-forge/main/control/src/sl_td04ad.cc
trunk/octave-forge/main/control/src/sl_tg01ad.cc
trunk/octave-forge/main/control/src/sl_tg01hd.cc
trunk/octave-forge/main/control/src/sl_tg01id.cc
trunk/octave-forge/main/control/src/sl_tg01jd.cc
trunk/octave-forge/main/control/src/sl_tg04bx.cc
Removed Paths:
-------------
trunk/octave-forge/main/control/src/slab01od.cc
trunk/octave-forge/main/control/src/slab04md.cc
trunk/octave-forge/main/control/src/slab08nd.cc
trunk/octave-forge/main/control/src/slab09hd.cc
trunk/octave-forge/main/control/src/slab09id.cc
trunk/octave-forge/main/control/src/slab09jd.cc
trunk/octave-forge/main/control/src/slab13ad.cc
trunk/octave-forge/main/control/src/slab13bd.cc
trunk/octave-forge/main/control/src/slab13dd.cc
trunk/octave-forge/main/control/src/slag08bd.cc
trunk/octave-forge/main/control/src/slib01ad.cc
trunk/octave-forge/main/control/src/slib01cd.cc
trunk/octave-forge/main/control/src/slident.cc
trunk/octave-forge/main/control/src/slsb01bd.cc
trunk/octave-forge/main/control/src/slsb02od.cc
trunk/octave-forge/main/control/src/slsb03md.cc
trunk/octave-forge/main/control/src/slsb03od.cc
trunk/octave-forge/main/control/src/slsb04md.cc
trunk/octave-forge/main/control/src/slsb04qd.cc
trunk/octave-forge/main/control/src/slsb10dd.cc
trunk/octave-forge/main/control/src/slsb10ed.cc
trunk/octave-forge/main/control/src/slsb10fd.cc
trunk/octave-forge/main/control/src/slsb10hd.cc
trunk/octave-forge/main/control/src/slsb10id.cc
trunk/octave-forge/main/control/src/slsb10jd.cc
trunk/octave-forge/main/control/src/slsb10kd.cc
trunk/octave-forge/main/control/src/slsb10yd.cc
trunk/octave-forge/main/control/src/slsb10zd.cc
trunk/octave-forge/main/control/src/slsb16ad.cc
trunk/octave-forge/main/control/src/slsb16bd.cc
trunk/octave-forge/main/control/src/slsb16cd.cc
trunk/octave-forge/main/control/src/slsg02ad.cc
trunk/octave-forge/main/control/src/slsg03ad.cc
trunk/octave-forge/main/control/src/slsg03bd.cc
trunk/octave-forge/main/control/src/sltb01id.cc
trunk/octave-forge/main/control/src/sltb01pd.cc
trunk/octave-forge/main/control/src/sltb01ud.cc
trunk/octave-forge/main/control/src/sltb04bd.cc
trunk/octave-forge/main/control/src/sltd04ad.cc
trunk/octave-forge/main/control/src/sltg01ad.cc
trunk/octave-forge/main/control/src/sltg01hd.cc
trunk/octave-forge/main/control/src/sltg01id.cc
trunk/octave-forge/main/control/src/sltg01jd.cc
trunk/octave-forge/main/control/src/sltg04bx.cc
Modified: trunk/octave-forge/main/control/src/__control_slicot_functions__.cc
===================================================================
--- trunk/octave-forge/main/control/src/__control_slicot_functions__.cc 2012-09-08 20:02:40 UTC (rev 10987)
+++ trunk/octave-forge/main/control/src/__control_slicot_functions__.cc 2012-09-09 19:09:48 UTC (rev 10988)
@@ -1,47 +1,47 @@
-#include "slab08nd.cc" // transmission zeros of state-space models
-#include "slab13dd.cc" // L-infinity norm
-#include "slsb10hd.cc" // H-2 controller synthesis - continuous-time
-#include "slsb10ed.cc" // H-2 controller synthesis - discrete-time
-#include "slab13bd.cc" // H-2 norm
-#include "slsb01bd.cc" // Pole assignment
-#include "slsb10fd.cc" // H-infinity controller synthesis - continuous-time
-#include "slsb10dd.cc" // H-infinity controller synthesis - discrete-time
-#include "slsb03md.cc" // Lyapunov equations
-#include "slsb04md.cc" // Sylvester equations - continuous-time
-#include "slsb04qd.cc" // Sylvester equations - discrete-time
-#include "slsg03ad.cc" // generalized Lyapunov equations
-#include "slsb02od.cc" // algebraic Riccati equations
-#include "slab13ad.cc" // Hankel singular values
-#include "slab01od.cc" // staircase form using orthogonal transformations
-#include "sltb01pd.cc" // minimal realization of state-space models
-#include "slsb03od.cc" // Cholesky factor of Lyapunov equations
-#include "slsg03bd.cc" // Cholesky factor of generalized Lyapunov equations
-#include "slag08bd.cc" // transmission zeros of descriptor state-space models
-#include "sltg01jd.cc" // minimal realization of descriptor state-space models
-#include "sltg01hd.cc" // controllability staircase form of descriptor state-space models
-#include "sltg01id.cc" // observability staircase form of descriptor state-space models
-#include "slsg02ad.cc" // solution of algebraic Riccati equations for descriptor systems
-#include "sltg04bx.cc" // gain of descriptor state-space models
-#include "sltb01id.cc" // scaling of state-space models
-#include "sltg01ad.cc" // scaling of descriptor state-space models
-#include "slsb10id.cc" // H-infinity loop shaping - continuous-time
-#include "slsb10kd.cc" // H-infinity loop shaping - discrete-time - strictly proper case
-#include "slsb10zd.cc" // H-infinity loop shaping - discrete-time - proper case
-#include "sltb04bd.cc" // State-space to transfer function conversion
-#include "slab04md.cc" // bilinear transformation
-#include "slsb10jd.cc" // descriptor to regular state-space conversion
-#include "sltd04ad.cc" // transfer function to state-space conversion
-#include "sltb01ud.cc" // controllable block Hessenberg realization
-#include "slab09hd.cc" // balanced stochastic truncation model reduction
-#include "slab09id.cc" // balanced truncation & singular perturbation approximation model reduction
-#include "slab09jd.cc" // hankel-norm approximation model reduction
-#include "slsb16ad.cc" // balanced truncation & singular perturbation approximation controller reduction
-#include "slsb16bd.cc" // coprime factorization state-feedback controller reduction
-#include "slsb16cd.cc" // frequency-weighted coprime factorization state-feedback controller reduction
-#include "slsb10yd.cc" // fit state-space model to frequency response data
-#include "slident.cc" // system identification
-#include "slib01cd.cc" // compute initial state vector
-#include "slib01ad.cc" // compute singular values
+#include "sl_ab08nd.cc" // transmission zeros of state-space models
+#include "sl_ab13dd.cc" // L-infinity norm
+#include "sl_sb10hd.cc" // H-2 controller synthesis - continuous-time
+#include "sl_sb10ed.cc" // H-2 controller synthesis - discrete-time
+#include "sl_ab13bd.cc" // H-2 norm
+#include "sl_sb01bd.cc" // Pole assignment
+#include "sl_sb10fd.cc" // H-infinity controller synthesis - continuous-time
+#include "sl_sb10dd.cc" // H-infinity controller synthesis - discrete-time
+#include "sl_sb03md.cc" // Lyapunov equations
+#include "sl_sb04md.cc" // Sylvester equations - continuous-time
+#include "sl_sb04qd.cc" // Sylvester equations - discrete-time
+#include "sl_sg03ad.cc" // generalized Lyapunov equations
+#include "sl_sb02od.cc" // algebraic Riccati equations
+#include "sl_ab13ad.cc" // Hankel singular values
+#include "sl_ab01od.cc" // staircase form using orthogonal transformations
+#include "sl_tb01pd.cc" // minimal realization of state-space models
+#include "sl_sb03od.cc" // Cholesky factor of Lyapunov equations
+#include "sl_sg03bd.cc" // Cholesky factor of generalized Lyapunov equations
+#include "sl_ag08bd.cc" // transmission zeros of descriptor state-space models
+#include "sl_tg01jd.cc" // minimal realization of descriptor state-space models
+#include "sl_tg01hd.cc" // controllability staircase form of descriptor state-space models
+#include "sl_tg01id.cc" // observability staircase form of descriptor state-space models
+#include "sl_sg02ad.cc" // solution of algebraic Riccati equations for descriptor systems
+#include "sl_tg04bx.cc" // gain of descriptor state-space models
+#include "sl_tb01id.cc" // scaling of state-space models
+#include "sl_tg01ad.cc" // scaling of descriptor state-space models
+#include "sl_sb10id.cc" // H-infinity loop shaping - continuous-time
+#include "sl_sb10kd.cc" // H-infinity loop shaping - discrete-time - strictly proper case
+#include "sl_sb10zd.cc" // H-infinity loop shaping - discrete-time - proper case
+#include "sl_tb04bd.cc" // State-space to transfer function conversion
+#include "sl_ab04md.cc" // bilinear transformation
+#include "sl_sb10jd.cc" // descriptor to regular state-space conversion
+#include "sl_td04ad.cc" // transfer function to state-space conversion
+#include "sl_tb01ud.cc" // controllable block Hessenberg realization
+#include "sl_ab09hd.cc" // balanced stochastic truncation model reduction
+#include "sl_ab09id.cc" // balanced truncation & singular perturbation approximation model reduction
+#include "sl_ab09jd.cc" // hankel-norm approximation model reduction
+#include "sl_sb16ad.cc" // balanced truncation & singular perturbation approximation controller reduction
+#include "sl_sb16bd.cc" // coprime factorization state-feedback controller reduction
+#include "sl_sb16cd.cc" // frequency-weighted coprime factorization state-feedback controller reduction
+#include "sl_sb10yd.cc" // fit state-space model to frequency response data
+#include "sl_ident.cc" // system identification
+#include "sl_ib01cd.cc" // compute initial state vector
+#include "sl_ib01ad.cc" // compute singular values
// stub function to avoid gen_doc_cache warning upon package installation
Copied: trunk/octave-forge/main/control/src/sl_ab01od.cc (from rev 10987, trunk/octave-forge/main/control/src/slab01od.cc)
===================================================================
--- trunk/octave-forge/main/control/src/sl_ab01od.cc (rev 0)
+++ trunk/octave-forge/main/control/src/sl_ab01od.cc 2012-09-09 19:09:48 UTC (rev 10988)
@@ -0,0 +1,139 @@
+/*
+
+Copyright (C) 2010 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 <http://www.gnu.org/licenses/>.
+
+Staircase controllability form.
+Uses SLICOT AB01OD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: August 2010
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ab01od, AB01OD)
+ (char& STAGES,
+ char& JOBU, char& JOBV,
+ int& N, int& M,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* U, int& LDU,
+ double* V, int& LDV,
+ int& NCONT, int& INDCON,
+ int* KSTAIR,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& INFO);
+}
+
+// PKG_ADD: autoload ("__sl_ab01od__", "__control_slicot_functions__.oct");
+DEFUN_DLD (__sl_ab01od__, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot AB01OD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 3)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char stages = 'F';
+ char jobu = 'I';
+ char jobv = 'N'; // not referenced because stages = 'F'
+
+ Matrix a = args(0).matrix_value ();
+ Matrix b = args(1).matrix_value ();
+ double tol = args(2).double_value ();
+
+ int n = a.rows (); // n: number of states
+ int m = b.columns (); // m: number of inputs
+
+ int lda = max (1, n);
+ int ldb = max (1, n);
+ int ldu = max (1, n);
+ int ldv = 1;
+
+ // arguments out
+ Matrix u (ldu, n);
+ double* v = 0; // not referenced because stages = 'F'
+
+ int ncont;
+ int indcon;
+
+ OCTAVE_LOCAL_BUFFER (int, kstair, n);
+
+ // workspace
+ int ldwork = max (1, n + max (n, 3*m));
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, m);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicators
+ int info;
+
+
+ // SLICOT routine AB01OD
+ F77_XFCN (ab01od, AB01OD,
+ (stages,
+ jobu, jobv,
+ n, m,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ u.fortran_vec (), ldu,
+ v, ldv,
+ ncont, indcon,
+ kstair,
+ tol,
+ iwork,
+ dwork, ldwork,
+ info));
+
+ if (f77_exception_encountered)
+ error ("__sl_ab01od__: exception in SLICOT subroutine AB01OD");
+
+ if (info != 0)
+ error ("__sl_ab01od__: AB01OD returned info = %d", info);
+
+ // resize
+ a.resize (n, n);
+ b.resize (n, m);
+ u.resize (n, n);
+
+ // return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = u;
+ retval(3) = octave_value (ncont);
+ }
+
+ return retval;
+}
Copied: trunk/octave-forge/main/control/src/sl_ab04md.cc (from rev 10987, trunk/octave-forge/main/control/src/slab04md.cc)
===================================================================
--- trunk/octave-forge/main/control/src/sl_ab04md.cc (rev 0)
+++ trunk/octave-forge/main/control/src/sl_ab04md.cc 2012-09-09 19:09:48 UTC (rev 10988)
@@ -0,0 +1,129 @@
+/*
+
+Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
+
+Discrete-time <--> continuous-time systems conversion
+by a bilinear transformation.
+Uses SLICOT AB04MD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: September 2011
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ab04md, AB04MD)
+ (char& TYPE,
+ int& N, int& M, int& P,
+ double& ALPHA, double& BETA,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* D, int& LDD,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& INFO);
+}
+
+// PKG_ADD: autoload ("__sl_ab04md__", "__control_slicot_functions__.oct");
+DEFUN_DLD (__sl_ab04md__, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot AB04MD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 7)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char type;
+
+ Matrix a = args(0).matrix_value ();
+ Matrix b = args(1).matrix_value ();
+ Matrix c = args(2).matrix_value ();
+ Matrix d = args(3).matrix_value ();
+
+ double alpha = args(4).double_value ();
+ double beta = args(5).double_value ();
+ int discrete = args(6).int_value ();
+
+ if (discrete == 0)
+ type = 'C';
+ else
+ type = 'D';
+
+ int n = a.rows (); // n: number of states
+ int m = b.columns (); // m: number of inputs
+ int p = c.rows (); // p: number of outputs
+
+ int lda = max (1, n);
+ int ldb = max (1, n);
+ int ldc = max (1, p);
+ int ldd = max (1, p);
+
+ // workspace
+ int ldwork = max (1, n);
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, n);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicator
+ int info;
+
+
+ // SLICOT routine AB04MD
+ F77_XFCN (ab04md, AB04MD,
+ (type,
+ n, m, p,
+ alpha, beta,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ c.fortran_vec (), ldc,
+ d.fortran_vec (), ldd,
+ iwork,
+ dwork, ldwork,
+ info));
+
+ if (f77_exception_encountered)
+ error ("__sl_ab04md__: exception in SLICOT subroutine AB04MD");
+
+ if (info != 0)
+ error ("__sl_ab04md__: AB04MD returned info = %d", info);
+
+ // return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = c;
+ retval(3) = d;
+ }
+
+ return retval;
+}
Copied: trunk/octave-forge/main/control/src/sl_ab08nd.cc (from rev 10987, trunk/octave-forge/main/control/src/slab08nd.cc)
===================================================================
--- trunk/octave-forge/main/control/src/sl_ab08nd.cc (rev 0)
+++ trunk/octave-forge/main/control/src/sl_ab08nd.cc 2012-09-09 19:09:48 UTC (rev 10988)
@@ -0,0 +1,225 @@
+/*
+
+Copyright (C) 2009, 2010, 2011 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 <http://www.gnu.org/licenses/>.
+
+Transmission zeros of state-space models.
+Uses SLICOT AB08ND by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: November 2009
+Version: 0.5
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+#include <complex>
+#include <xpow.h>
+
+extern "C"
+{
+ int F77_FUNC (ab08nd, AB08ND)
+ (char& EQUIL,
+ int& N, int& M, int& P,
+ const double* A, int& LDA,
+ const double* B, int& LDB,
+ const double* C, int& LDC,
+ const double* D, int& LDD,
+ int& NU, int& RANK, int& DINFZ,
+ int& NKROR, int& NKROL, int* INFZ,
+ int* KRONR, int* KRONL,
+ double* AF, int& LDAF,
+ double* BF, int& LDBF,
+ double& TOL,
+ int* IWORK, double* DWORK, int& LDWORK,
+ int& INFO);
+
+ int F77_FUNC (dggev, DGGEV)
+ (char& JOBVL, char& JOBVR,
+ int& N,
+ double* AF, int& LDAF,
+ double* BF, int& LDBF,
+ double* ALPHAR, double* ALPHAI,
+ double* BETA,
+ double* VL, int& LDVL,
+ double* VR, int& LDVR,
+ double* WORK, int& LWORK,
+ int& INFO);
+}
+
+// PKG_ADD: autoload ("__sl_ab08nd__", "__control_slicot_functions__.oct");
+DEFUN_DLD (__sl_ab08nd__, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot AB08ND Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 5)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char equil;
+
+ const Matrix a = args(0).matrix_value ();
+ const Matrix b = args(1).matrix_value ();
+ const Matrix c = args(2).matrix_value ();
+ const Matrix d = args(3).matrix_value ();
+ const int scaled = args(4).int_value ();
+
+ if (scaled == 0)
+ equil = 'S';
+ else
+ equil = 'N';
+
+ int n = a.rows (); // n: number of states
+ int m = b.columns (); // m: number of inputs
+ int p = c.rows (); // p: number of outputs
+
+ int lda = max (1, a.rows ());
+ int ldb = max (1, b.rows ());
+ int ldc = max (1, c.rows ());
+ int ldd = max (1, d.rows ());
+
+ // arguments out
+ int nu;
+ int rank;
+ int dinfz;
+ int nkror;
+ int nkrol;
+
+ int ldaf = max (1, n + m);
+ int ldbf = max (1, n + p);
+
+ OCTAVE_LOCAL_BUFFER (int, infz, n);
+ OCTAVE_LOCAL_BUFFER (int, kronr, 1 + max (n, m));
+ OCTAVE_LOCAL_BUFFER (int, kronl, 1 + max (n, p));
+
+ OCTAVE_LOCAL_BUFFER (double, af, ldaf * (n + min (p, m)));
+ OCTAVE_LOCAL_BUFFER (double, bf, ldbf * (n + m));
+
+ // workspace
+ int s = max (m, p);
+ int ldwork = max (s, n) + max (3*s-1, n+s);
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, s);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicator
+ int info;
+
+ // tolerance
+ double tol = 0; // AB08ND uses DLAMCH for default tolerance
+
+ // SLICOT routine AB08ND
+ F77_XFCN (ab08nd, AB08ND,
+ (equil,
+ n, m, p,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ c.fortran_vec (), ldc,
+ d.fortran_vec (), ldd,
+ nu, rank, dinfz,
+ nkror, nkrol, infz,
+ kronr, kronl,
+ af, ldaf,
+ bf, ldbf,
+ tol,
+ iwork, dwork, ldwork,
+ info));
+
+ if (f77_exception_encountered)
+ error ("ss: zero: __sl_ab08nd__: exception in SLICOT subroutine AB08ND");
+
+ if (info != 0)
+ error ("ss: zero: __sl_ab08nd__: AB08ND returned info = %d", info);
+
+
+ // DGGEV Part
+ char jobvl = 'N';
+ char jobvr = 'N';
+
+ double* vl = 0; // not referenced because jobvl = 'N'
+ int ldvl = 1;
+ double* vr = 0; // not referenced because jobvr = 'N'
+ int ldvr = 1;
+
+ int lwork = max (1, 8*nu);
+ OCTAVE_LOCAL_BUFFER (double, work, lwork);
+
+ ColumnVector alphar (nu);
+ ColumnVector alphai (nu);
+ ColumnVector beta (nu);
+
+ int info2;
+
+ F77_XFCN (dggev, DGGEV,
+ (jobvl, jobvr,
+ nu,
+ af, ldaf,
+ bf, ldbf,
+ alphar.fortran_vec (), alphai.fortran_vec (),
+ beta.fortran_vec (),
+ vl, ldvl,
+ vr, ldvr,
+ work, lwork,
+ info2));
+
+ if (f77_exception_encountered)
+ error ("ss: zero: __sl_ab08nd__: exception in LAPACK subroutine DGGEV");
+
+ if (info2 != 0)
+ error ("ss: zero: __sl_ab08nd__: DGGEV returned info = %d", info2);
+
+ // calculate gain
+ octave_value gain = Matrix (0, 0);;
+
+ if (m == 1 && p == 1)
+ {
+ if (nu < n)
+ gain = c * xpow (a, double (n-1-nu)) * b;
+ else
+ gain = d;
+ }
+
+ // assemble complex vector - adapted from DEFUN complex in data.cc
+ ColumnVector zeror (nu);
+ ColumnVector zeroi (nu);
+
+ zeror = quotient (alphar, beta);
+ zeroi = quotient (alphai, beta);
+
+ ComplexColumnVector zero (nu, Complex ());
+
+ for (octave_idx_type i = 0; i < nu; i++)
+ zero.xelem (i) = Complex (zeror(i), zeroi(i));
+
+ // return values
+ retval(0) = zero;
+ retval(1) = gain;
+ }
+
+ return retval;
+}
Copied: trunk/octave-forge/main/control/src/sl_ab09hd.cc (from rev 10987, trunk/octave-forge/main/control/src/slab09hd.cc)
===================================================================
--- trunk/octave-forge/main/control/src/sl_ab09hd.cc (rev 0)
+++ trunk/octave-forge/main/control/src/sl_ab09hd.cc 2012-09-09 19:09:48 UTC (rev 10988)
@@ -0,0 +1,242 @@
+/*
+
+Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
+
+Model reduction based on balanced stochastic truncation method.
+Uses SLICOT AB09HD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: October 2011
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ab09hd, AB09HD)
+ (char& DICO, char& JOB, char& EQUIL, char& ORDSEL,
+ int& N, int& M, int& P,
+ int& NR,
+ double& ALPHA, double& BETA,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* D, int& LDD,
+ int& NS,
+ double* HSV,
+ double& TOL1, double& TOL2,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ bool* BWORK,
+ int& IWARN, int& INFO);
+}
+
+// PKG_ADD: autoload ("__sl_ab09hd__", "__control_slicot_functions__.oct");
+DEFUN_DLD (__sl_ab09hd__, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot AB09HD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 13)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char dico;
+ char job;
+ char equil;
+ char ordsel;
+
+ Matrix a = args(0).matrix_value ();
+ Matrix b = args(1).matrix_value ();
+ Matrix c = args(2).matrix_value ();
+ Matrix d = args(3).matrix_value ();
+
+ const int idico = args(4).int_value ();
+ const int iequil = args(5).int_value ();
+ const int ijob = args(6).int_value ();
+
+ int nr = args(7).int_value ();
+ const int iordsel = args(8).int_value ();
+
+ double alpha = args(9).double_value ();
+ double beta = args(10).double_value ();
+ double tol1 = args(11).double_value ();
+ double tol2 = args(12).double_value ();
+
+ switch (ijob)
+ {
+ case 0:
+ job = 'B';
+ break;
+ case 1:
+ job = 'F';
+ break;
+ case 2:
+ job = 'S';
+ break;
+ case 3:
+ job = 'P';
+ break;
+ default:
+ error ("__sl_ab09hd__: argument job invalid");
+ }
+
+ if (idico == 0)
+ dico = 'C';
+ else
+ dico = 'D';
+
+ if (iequil == 0)
+ equil = 'S';
+ else
+ equil = 'N';
+
+ if (iordsel == 0)
+ ordsel = 'F';
+ else
+ ordsel = 'A';
+
+ int n = a.rows (); // n: number of states
+ int m = b.columns (); // m: number of inputs
+ int p = c.rows (); // p: number of outputs
+
+ int lda = max (1, n);
+ int ldb = max (1, n);
+ int ldc = max (1, p);
+ int ldd = max (1, p);
+
+ // arguments out
+ int ns;
+ ColumnVector hsv (n);
+
+ // workspace
+ int liwork = max (1, 2*n);
+ int mb;
+
+ if (beta == 0)
+ mb = m;
+ else
+ mb = m + p;
+
+ int ldwork = 2*n*n + mb*(n+p)
+ + max (2, n*(max (n,mb,p)+5), 2*n*p + max (p*(mb+2), 10*n*(n+1)));
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+ OCTAVE_LOCAL_BUFFER (bool, bwork, 2*n);
+
+ // error indicators
+ int iwarn = 0;
+ int info = 0;
+
+
+ // SLICOT routine AB09HD
+ F77_XFCN (ab09hd, AB09HD,
+ (dico, job, equil, ordsel,
+ n, m, p,
+ nr,
+ alpha, beta,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ c.fortran_vec (), ldc,
+ d.fortran_vec (), ldd,
+ ns,
+ hsv.fortran_vec (),
+ tol1, tol2,
+ iwork,
+ dwork, ldwork,
+ bwork,
+ iwarn, info));
+
+ if (f77_exception_encountered)
+ error ("bstmodred: exception in SLICOT subroutine AB09HD");
+
+
+ static const char* err_msg[] = {
+ "0: OK",
+ "1: the computation of the ordered real Schur form of A "
+ "failed",
+ "2: the reduction of the Hamiltonian matrix to real "
+ "Schur form failed",
+ "3: the reordering of the real Schur form of the "
+ "Hamiltonian matrix failed",
+ "4: the Hamiltonian matrix has less than N stable "
+ "eigenvalues",
+ "5: the coefficient matrix U11 in the linear system "
+ "X*U11 = U21 to determine X is singular to working "
+ "precision",
+ "6: BETA = 0 and D has not a maximal row rank",
+ "7: the computation of Hankel singular values failed",
+ "8: the separation of the ALPHA-stable/unstable diagonal "
+ "blocks failed because of very close eigenvalues",
+ "9: the resulting order of reduced stable part is less "
+ "than the number of unstable zeros of the stable "
+ "part"};
+
+ static const char* warn_msg[] = {
+ "0: OK",
+ "1: with ORDSEL = 'F', the selected order NR is greater "
+ "than NSMIN, the sum of the order of the "
+ "ALPHA-unstable part and the order of a minimal "
+ "realization of the ALPHA-stable part of the given "
+ "system; in this case, the resulting NR is set equal "
+ "to NSMIN.",
+ "2: with ORDSEL = 'F', the selected order NR corresponds "
+ "to repeated singular values for the ALPHA-stable "
+ "part, which are neither all included nor all "
+ "excluded from the reduced model; in this case, the "
+ "resulting NR is automatically decreased to exclude "
+ "all repeated singular values.",
+ "3: with ORDSEL = 'F', the selected order NR is less "
+ "than the order of the ALPHA-unstable part of the "
+ "given system; in this case NR is set equal to the "
+ "order of the ALPHA-unstable part."};
+
+ error_msg ("bstmodred", info, 9, err_msg);
+ warning_msg ("bstmodred", iwarn, 3, warn_msg);
+
+ // resize
+ a.resize (nr, nr);
+ b.resize (nr, m);
+ c.resize (p, nr);
+ hsv.resize (ns);
+
+ // return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = c;
+ retval(3) = d;
+ retval(4) = octave_value (nr);
+ retval(5) = hsv;
+ retval(6) = octave_value (ns);
+ }
+
+ return retval;
+}
Copied: trunk/octave-forge/main/control/src/sl_ab09id.cc (from rev 10987, trunk/octave-forge/main/control/src/slab09id.cc)
===================================================================
--- trunk/octave-forge/main/control/src/sl_ab09id.cc (rev 0)
+++ trunk/octave-forge/main/control/src/sl_ab09id.cc 2012-09-09 19:09:48 UTC (rev 10988)
@@ -0,0 +1,396 @@
+/*
+
+Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
+
+Model reduction based on Balance & Truncate (B&T) or
+Singular Perturbation Approximation (SPA) method.
+Uses SLICOT AB09ID by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: October 2011
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ab09id, AB09ID)
+ (char& DICO, char& JOBC, char& JOBO, char& JOB,
+ char& WEIGHT, char& EQUIL, char& ORDSEL,
+ int& N, int& M, int& P,
+ int& NV, int& PV, int& NW, int& MW,
+ int& NR,
+ double& ALPHA, double& ALPHAC, double& ALPHAO,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* D, int& LDD,
+ double* AV, int& LDAV,
+ double* BV, int& LDBV,
+ double* CV, int& LDCV,
+ double* DV, int& LDDV,
+ double* AW, int& LDAW,
+ double* BW, int& LDBW,
+ double* CW, int& LDCW,
+ double* DW, int& LDDW,
+ int& NS,
+ double* HSV,
+ double& TOL1, double& TOL2,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+}
+
+// PKG_ADD: autoload ("__sl_ab09id__", "__control_slicot_functions__.oct");
+DEFUN_DLD (__sl_ab09id__, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot AB09ID Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 25)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char dico;
+ char jobc;
+ char jobo;
+ char job;
+ char weight;
+ char equil;
+ char ordsel;
+
+ Matrix a = args(0).matrix_value ();
+ Matrix b = args(1).matrix_value ();
+ Matrix c = args(2).matrix_value ();
+ Matrix d = args(3).matrix_value ();
+
+ const int idico = args(4).int_value ();
+ const int iequil = args(5).int_value ();
+ int nr = args(6).int_value ();
+ const int iordsel = args(7).int_value ();
+ double alpha = args(8).double_value ();
+ const int ijob = args(9).int_value ();
+
+ Matrix av = args(10).matrix_value ();
+ Matrix bv = args(11).matrix_value ();
+ Matrix cv = args(12).matrix_value ();
+ Matrix dv = args(13).matrix_value ();
+
+ Matrix aw = args(14).matrix_value ();
+ Matrix bw = args(15).matrix_value ();
+ Matrix cw = args(16).matrix_value ();
+ Matrix dw = args(17).matrix_value ();
+
+ const int iweight = args(18).int_value ();
+ const int ijobc = args(19).int_value ();
+ double alphac = args(20).double_value ();
+ const int ijobo = args(21).int_value ();
+ double alphao = args(22).double_value ();
+
+ double tol1 = args(23).double_value ();
+ double tol2 = args(24).double_value ();
+
+ if (idico == 0)
+ dico = 'C';
+ else
+ dico = 'D';
+
+ if (iequil == 0)
+ equil = 'S';
+ else
+ equil = 'N';
+
+ if (iordsel == 0)
+ ordsel = 'F';
+ else
+ ordsel = 'A';
+
+ if (ijobc == 0)
+ jobc = 'S';
+ else
+ jobc = 'E';
+
+ if (ijobo == 0)
+ jobo = 'S';
+ else
+ jobo = 'E';
+
+ switch (ijob)
+ {
+ case 0:
+ job = 'B';
+ break;
+ case 1:
+ job = 'F';
+ break;
+ case 2:
+ job = 'S';
+ break;
+ case 3:
+ job = 'P';
+ break;
+ default:
+ error ("__sl_ab09id__: argument job invalid");
+ }
+
+ switch (iweight)
+ {
+ case 0:
+ weight = 'N';
+ break;
+ case 1:
+ weight = 'L';
+ break;
+ case 2:
+ weight = 'R';
+ break;
+ case 3:
+ weight = 'B';
+ break;
+ default:
+ error ("__sl_ab09id__: argument weight invalid");
+ }
+
+ int n = a.rows (); // n: number of states
+ int m = b.columns (); // m: number of inputs
+ int p = c.rows (); // p: number of outputs
+
+ int nv = av.rows ();
+ int pv = cv.rows ();
+ int nw = aw.rows ();
+ int mw = bw.columns ();
+
+ int lda = max (1, n);
+ int ldb = max (1, n);
+ int ldc = max (1, p);
+ int ldd = max (1, p);
+
+ int ldav = max (1, nv);
+ int ldbv = max (1, nv);
+ int ldcv = max (1, pv);
+ int lddv = max (1, pv);
+
+ int ldaw = max (1, nw);
+ int ldbw = max (1, nw);
+ int ldcw = max (1, m);
+ int lddw = max (1, m);
+
+ // arguments out
+ int ns;
+ ColumnVector hsv (n);
+
+ // workspace
+ int liwork;
+ int liwrk1;
+ int liwrk2;
+ int liwrk3;
+
+ switch (job)
+ {
+ case 'B':
+ liwrk1 = 0;
+ break;
+ case 'F':
+ liwrk1 = n;
+ break;
+ default:
+ liwrk1 = 2*n;
+ }
+
+ if (nv == 0 || weight == 'R' || weight == 'N')
+ liwrk2 = 0;
+ else
+ liwrk2 = nv + max (p, pv);
+
+ if (nw == 0 || weight == 'L' || weight == 'N')
+ liwrk3 = 0;
+ else
+ liwrk3 = nw + max (m, mw);
+
+ liwork = max (3, liwrk1, liwrk2, liwrk3);
+
+ int ldwork;
+ int lminl;
+ int lrcf;
+ int lminr;
+ int llcf;
+ int lleft;
+ int lright;
+
+ if (nw == 0 || weight == 'L' || weight == 'N')
+ {
+ lrcf = 0;
+ lminr = 0;
+ }
+ else
+ {
+ lrcf = mw*(nw+mw) + max (nw*(nw+5), mw*(mw+2), 4*mw, 4*m);
+ if (m == mw)
+ lminr = nw + max (nw, 3*m);
+ else
+ lminr = 2*nw*max (m, mw) + nw + max (nw, 3*m, 3*mw);
+ }
+
+ llcf = pv*(nv+pv) + pv*nv + max (nv*(nv+5), pv*(pv+2), 4*pv, 4*p);
+
+ if (nv == 0 || weight == 'R' || weight == 'N')
+ lminl = 0;
+ else if (p == pv)
+ lminl = max (llcf, nv + max (nv, 3*p));
+ else
+ lminl = max (p, pv) * (2*nv + max (p, pv)) + max (llcf, nv + max (nv, 3*p, 3*pv));
+
+
+ if (pv == 0 || weight == 'R' || weight == 'N')
+ lleft = n*(p+5);
+ else
+ lleft = (n+nv) * (n + nv + max (n+nv, pv) + 5);
+
+ if (mw == 0 || weight == 'L' || weight == 'N')
+ lright = n*(m+5);
+ else
+ lright = (n+nw) * (n + nw + max (n+nw, mw) + 5);
+
+ ldwork = max (lminl, lminr, lrcf,
+ 2*n*n + max (1, lleft, lright, 2*n*n+5*n, n*max (m, p)));
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicators
+ int iwarn = 0;
+ int info = 0;
+
+
+ // SLICOT routine AB09ID
+ F77_XFCN (ab09id, AB09ID,
+ (dico, jobc, jobo, job,
+ weight, equil, ordsel,
+ n, m, p,
+ nv, pv, nw, mw,
+ nr,
+ alpha, alphac, alphao,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ c.fortran_vec (), ldc,
+ d.fortran_vec (), ldd,
+ av.fortran_vec (), ldav,
+ bv.fortran_vec (), ldbv,
+ cv.fortran_vec (), ldcv,
+ dv.fortran_vec (), lddv,
+ aw.fortran_vec (), ldaw,
+ bw.fortran_vec (), ldbw,
+ cw.fortran_vec (), ldcw,
+ dw.fortran_vec (), lddw,
+ ns,
+ hsv.fortran_vec (),
+ tol1, tol2,
+ iwork,
+ dwork, ldwork,
+ iwarn, info));
+
+ if (f77_exception_encountered)
+ error ("modred: exception in SLICOT subroutine AB09ID");
+
+
+ static const char* err_msg[] = {
+ "0: OK",
+ "1: the computation of the ordered real Schur form of A "
+ "failed",
+ "2: the separation of the ALPHA-stable/unstable "
+ "diagonal blocks failed because of very close "
+ "eigenvalues",
+ "3: the reduction to a real Schur form of the state "
+ "matrix of a minimal realization of V failed",
+ "4: a failure was detected during the ordering of the "
+ "real Schur form of the state matrix of a minimal "
+ "realization of V or in the iterative process to "
+ "compute a left coprime factorization with inner "
+ "denominator",
+ "5: if DICO = 'C' and the matrix AV has an observable "
+ "eigenvalue on the imaginary axis, or DICO = 'D' and "
+ "AV has an observable eigenvalue on the unit circle",
+ "6: the reduction to a real Schur form of the state "
+ "matrix of a minimal realization of W failed",
+ "7: a failure was detected during the ordering of the "
+ "real Schur form of the state matrix of a minimal "
+ "realization of W or in the iterative process to "
+ "compute a right coprime factorization with inner "
+ "denominator",
+ "8: if DICO = 'C' and the matrix AW has a controllable "
+ "eigenvalue on the imaginary axis, or DICO = 'D' and "
+ "AW has a controllable eigenvalue on the unit circle",
+ "9: the computation of eigenvalues failed",
+ "10: the computation of Hankel singular values failed"};
+
+ static const char* warn_msg[] = {
+ "0: OK",
+ "1: with ORDSEL = 'F', the selected order NR is greater "
+ "than NSMIN, the sum of the order of the "
+ "ALPHA-unstable part and the order of a minimal "
+ "realization of the ALPHA-stable part of the given "
+ "system; in this case, the resulting NR is set equal "
+ "to NSMIN.",
+ "2: with ORDSEL = 'F', the selected order NR corresponds "
+ "to repeated singular values for the ALPHA-stable "
+ "part, which are neither all included nor all "
+ "excluded from the reduced model; in this case, the "
+ "resulting NR is automatically decreased to exclude "
+ "all repeated singular values.",
+ "3: with ORDSEL = 'F', the selected order NR is less "
+ "than the order of the ALPHA-unstable part of the "
+ "given system; in this case NR is set equal to the "
+ "order of the ALPHA-unstable part.",
+/* 10+%d: %d */ "violations of the numerical stability condition "
+ "occured during the assignment of eigenvalues in the "
+ "SLICOT Library routines SB08CD and/or SB08DD."};
+
+
+ error_msg ("modred", info, 10, err_msg);
+ warning_msg ("modred", iwarn, 3, warn_msg, 10);
+
+ // resize
+ a.resize (nr, nr);
+ b.resize (nr, m);
+ c.resize (p, nr);
+ hsv.resize (ns);
+
+ // return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = c;
+ retval(3) = d;
+ retval(4) = octave_value (nr);
+ retval(5) = hsv;
+ retval(6) = octave_value (ns);
+ }
+
+ return retval;
+}
Copied: trunk/octave-forge/main/control/src/sl_ab09jd.cc (from rev 10987, trunk/octave-forge/main/control/src/slab09jd.cc)
===================================================================
--- trunk/octave-forge/main/control/src/sl_ab09jd.cc (rev 0)
+++ trunk/octave-forge/main/control/src/sl_ab09jd.cc 2012-09-09 19:09:48 UTC (rev 10988)
@@ -0,0 +1,404 @@
+/*
+
+Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
+
+Model reduction based on Hankel-norm approximation method.
+Uses SLICOT AB09JD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: July 2011
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ab09jd, AB09JD)
+ (char& JOBV, char& JOBW, char& JOBINV,
+ char& DICO, char& EQUIL, char& ORDSEL,
+ int& N, int& NV, int& NW, int& M, int& P,
+ int& NR,
+ double& ALPHA,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* D, int& LDD,
+ double* AV, int& LDAV,
+ double* BV, int& LDBV,
+ double* CV, int& LDCV,
+ double* DV, int& LDDV,
+ double* AW, int& LDAW,
+ double* BW, int& LDBW,
+ double* CW, int& LDCW,
+ double* DW, int& LDDW,
+ int& NS,
+ double* HSV,
+ double& TOL1, double& TOL2,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+}
+
+// PKG_ADD: autoload ("__sl_ab09jd__", "__control_slicot_functions__.oct");
+DEFUN_DLD (__sl_ab09jd__, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot AB09JD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 22)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char jobv;
+ char jobw;
+ char jobinv;
+ char dico;
+ char equil;
+ char ordsel;
+
+ Matrix a = args(0).matrix_value ();
+ Matrix b = args(1).matrix_value ();
+ Matrix c = args(2).matrix_value ();
+ Matrix d = args(3).matrix_value ();
+
+ const int idico = args(4).int_value ();
+ const int iequil = args(5).int_value ();
+ int nr = args(6).int_value ();
+ const int iordsel = args(7).int_value ();
+ double alpha = args(8).double_value ();
+
+ const int ijobv = args(9).int_value ();
+ Matrix av = args(10).matrix_value ();
+ Matrix bv = args(11).matrix_value ();
+ Matrix cv = args(12).matrix_value ();
+ Matrix dv = args(13).matrix_value ();
+
+ const int ijobw = args(14).int_value ();
+ Matrix aw = args(15).matrix_value ();
+ Matrix bw = args(16).matrix_value ();
+ Matrix cw = args(17).matrix_value ();
+ Matrix dw = args(18).matrix_value ();
+
+ const int ijobinv = args(19).int_value ();
+ double tol1 = args(20).double_value ();
+ double tol2 = args(21).double_value ();
+
+ switch (ijobv)
+ {
+ case 0:
+ jobv = 'N';
+ break;
+ case 1:
+ jobv = 'V';
+ break;
+ case 2:
+ jobv = 'I';
+ break;
+ case 3:
+ jobv = 'C';
+ break;
+ case 4:
+ jobv = 'R';
+ break;
+ default:
+ error ("__sl_ab09jd__: argument jobv invalid");
+ }
+
+ switch (ijobw)
+ {
+ case 0:
+ jobw = 'N';
+ break;
+ case 1:
+ jobw = 'W';
+ break;
+ case 2:
+ jobw = 'I';
+ break;
+ case 3:
+ jobw = 'C';
+ break;
+ case 4:
+ jobw = 'R';
+ break;
+ default:
+ error ("__sl_ab09jd__: argument jobw invalid");
+ }
+
+ switch (ijobinv)
+ {
+ case 0:
+ jobinv = 'N';
+ break;
+ case 1:
+ jobinv = 'I';
+ break;
+ case 2:
+ jobinv = 'A';
+ break;
+ default:
+ error ("__sl_ab09jd__: argument jobinv invalid");
+ }
+
+ if (idico == 0)
+ dico = 'C';
+ else
+ dico = 'D';
+
+ if (iequil == 0)
+ equil = 'S';
+ else
+ equil = 'N';
+
+ if (iordsel == 0)
+ ordsel = 'F';
+ else
+ ordsel = 'A';
+
+ int n = a.rows (); // n: number of states
+ int nv = av.rows ();
+ int nw = aw.rows ();
+ int m = b.columns (); // m: number of inputs
+ int p = c.rows (); // p: number of outputs
+
+ int lda = max (1, n);
+ int ldb = max (1, n);
+ int ldc = max (1, p);
+ int ldd = max (1, p);
+
+ int ldav = max (1, nv);
+ int ldbv = max (1, nv);
+ int ldcv = max (1, p);
+ int lddv = max (1, p);
+
+ int ldaw = max (1, nw);
+ int ldbw = max (1, nw);
+ int ldcw = max (1, m);
+ int lddw = max (1, m);
+
+ // arguments out
+ int ns;
+ ColumnVector hsv (n);
+
+ // workspace
+ int liwork;
+ int tmpc;
+ int tmpd;
+
+ if (jobv == 'N')
+ tmpc = 0;
+ else
+ tmpc = max (2*p, nv+p+n+6, 2*nv+p+2);
+
+ if (jobw == 'N')
+ tmpd = 0;
+ else
+ tmpd = max (2*m, nw+m+n+6, 2*nw+m+2);
+
+ if (dico == 'C')
+ liwork = max (1, m, tmpc, tmpd);
+ else
+ liwork = max (1, n, m, tmpc, tmpd);
+
+ int ldwork;
+ int nvp = nv + p;
+ int nwm = nw + m;
+ int ldw1;
+ int ldw2;
+ int ldw3 = n*(2*n + max (n, m, p) + 5) + n*(n+1)/2;
+ int ldw4 = n*(m+p+2) + 2*m*p + min (n, m) + max (3*m+1, min (n, m) + p);
+
+ if (jobv == 'N')
+ {
+ ldw1 = 0;
+ }
+ else
+ {
+ ldw1 = 2*nvp*(nvp+p) + p*p + max (2*nvp*nvp + max (11*nvp+16, p*nvp),
+ nvp*n + max (nvp*n+n*n, p*n, p*m));
+ }
+
+ if (jobw == 'N')
+ {
+ ldw2 = 0;
+ }
+ else
+ {
+ ldw2 = 2*nwm*(nwm+m) + m*m + max (2*nwm*nwm + max (11*nwm+16, m*nwm),
+ nwm*n + max (nwm*n+n*n, m*n, p*m));
+ }
+
+ ldwork = max (ldw1, ldw2, ldw3, ldw4);
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicators
+ int iwarn = 0;
+ int info = 0;
+
+
+ // SLICOT routine AB09JD
+ F77_XFCN (ab09jd, AB09JD,
+ (jobv, jobw, jobinv,
+ dico, equil, ordsel,
+ n, nv, nw, m, p,
+ nr,
+ alpha,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ c.fortran_vec (), ldc,
+ d.fortran_vec (), ldd,
+ av.fortran_vec (), ldav,
+ bv.fortran_vec (), ldbv,
+ cv.fortran_vec (), ldcv,
+ dv.fortran_vec (), lddv,
+ aw.fortran_vec (), ldaw,
+ bw.fortran_vec (), ldbw,
+ cw.fortran_vec (), ldcw,
+ dw.fortran_vec (), lddw,
+ ns,
+ hsv.fortran_vec (),
+ tol1, tol2,
+ iwork,
+ dwork, ldwork,
+ iwarn, info));
+
+ if (f77_exception_encountered)
+ error ("hnamodred: exception in SLICOT subroutine AB09JD");
+
+
+ static const char* err_msg[] = {
+ "0: OK",
+ "1: the computation of the ordered real Schur form of A "
+ "failed",
+
+ "2: the separation of the ALPHA-stable/unstable "
+ "diagonal blocks failed because of very close eigenvalues",
+
+ "3: the reduction of AV to a real Schur form failed",
+
+ "4: the reduction of AW to a real Schur form failed",
+
+ "5: the reduction to generalized Schur form of the "
+ "descriptor pair corresponding to the inverse of V "
+ "failed",
+
+ "6: the reduction to generalized Schur form of the "
+ "descriptor pair corresponding to the inverse of W "
+ "failed",
+
+ "7: the computation of Hankel singular values failed",
+
+ "8: the computation of stable projection in the "
+ "Hankel-norm approximation algorithm failed",
+
+ "9: the order of computed stable projection in the "
+ "Hankel-norm approximation algorithm differs "
+ "from the order of Hankel-norm approximation",
+
+ "10: the reduction of AV-BV*inv(DV)*CV to a "
+ "real Schur form failed",
+
+ "11: the reduction of AW-BW*inv(DW)*CW to a "
+ "real Schur form failed",
+
+ "12: the solution of the Sylvester equation failed "
+ "because the poles of V (if JOBV = 'V') or of "
+ "conj(V) (if JOBV = 'C') are not distinct from "
+ "the poles of G1 (see METHOD)",
+
+ "13: the solution of the Sylvester equation failed "
+ "because the poles of W (if JOBW = 'W') or of "
+ "conj(W) (if JOBW = 'C') are not distinct from "
+ "the poles of G1 (see METHOD)",
+
+ "14: the solution of the Sylvester equation failed "
+ "because the zeros of V (if JOBV = 'I') or of "
+ "conj(V) (if JOBV = 'R') are not distinct from "
+ "the poles of G1sr (see METHOD)",
+
+ "15: the solution of the Sylvester equation failed "
+ "because the zeros of W (if JOBW = 'I') or of "
+ "conj(W) (if JOBW = 'R') are not distinct from "
+ "the poles of G1sr (see METHOD)",
+
+ "16: the solution of the generalized Sylvester system "
+ "failed because the zeros of V (if JOBV = 'I') or "
+ "of conj(V) (if JOBV = 'R') are not distinct from "
+ "the poles of G1sr (see METHOD)",
+
+ "17: the solution of the generalized Sylvester system "
+ "failed becau...
[truncated message content] |
