1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2// Copyright (C) INRIA -
3//
4// Copyright (C) 2012 - 2016 - Scilab Enterprises
5//
6// This file is hereby licensed under the terms of the GNU GPL v2.0,
7// pursuant to article 5.3.4 of the CeCILL v.2.1.
8// This file was originally licensed under the terms of the CeCILL v2.1,
9// and continues to be available under such terms.
10// For more information, see the COPYING file which you should have received
11// along with this program.
12
13function [Inn,X,Gbar]=rowinout(G)
14    // Inner-outer factorization (and row compression) of (lxp) G =:[A,B,C,D] with l>=p
15    // G is assumed to be tall (l>=p) without zero on the imaginary axis
16    // and with a D matrix which is full column rank. G must also be stable for
17    // having Gbar stable.
18    //
19    // G admits the following inner-outer factorization:
20    //
21    //          G = [ Inn ] | Gbar |
22    //                      |  0   |
23    //
24    // where Inn is square and inner (all pass and stable) and Gbar square and outer i.e:
25    // Gbar is square bi-proper and bi-stable (Gbar inverse is also proper and stable);
26    //
27    // Note that:
28    //          [ Gbar ]
29    //    X*G = [  -   ]
30    //          [  0   ]
31    // is a row compression of G where X = Inn inverse is all-pass:
32    //               T
33    // X is lxl and X (-s) X(s) = Identity (all-pass property).
34    //
35
36    G1=G(1);
37    flag="ss";if G1(1)=="r" then flag="tf";G=tf2ss(G);end
38    [rows,cols]=size(G);
39    [A,B,C,D]=abcd(G);
40    sqD=inv(sqrtm(D'*D));
41    //----------------
42    [w,rt]=rowcomp(D);
43    Dorto=w(rt+1:rows,:)';
44    // Inner factor:
45    //-------------
46    [F,Xc]=lqr(G);
47    Af=A+B*F;Cf=C+D*F;
48    Inn=syslin("c",Af,[B*sqD,-pinv(Xc)*C'*Dorto],Cf,[D*sqD,Dorto]);
49    X=invsyslin(Inn);
50    //----------------------------
51    // Outer factor:
52    Gbar=invsyslin(syslin("c",Af,B*sqD,F,sqD));
53    if flag=="tf" then Inner=ss2tf(Inner);Gbar=ss2tf(Gbar);X=ss2tf(X);end
54endfunction
55