1
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
4//
5// Copyright (C) 2012 - 2016 - Scilab Enterprises
6//
7// This file is hereby licensed under the terms of the GNU GPL v2.0,
8// pursuant to article 5.3.4 of the CeCILL v.2.1.
9// This file was originally licensed under the terms of the CeCILL v2.1,
10// and continues to be available under such terms.
11// For more information, see the COPYING file which you should have received
12// along with this program.
13
14function [Xp,dima,dimb,dim]=spantwo(A,B)
15    //Given two matrices A and B with same number of rows,
16    //returns a square matrix Xp (not necessarily orthogonal)
17    //such that :
18    //         [*,0]    (dim-dimb rows)
19    //Xp*[A,B]=[*,*]    (dima+dimb-dim rows)
20    //         [0,*]    (dim-dima rows)
21    //         [0,0]
22    //The  dima first columns of inv(Xp) span range(A).
23    //Columns dim-dimb+1 to dima of inv(Xp) span the intesection of
24    //range(A) and range(B).
25    //The dim first columns of inv(Xp) span range(A)+range(B).
26    // Ex: A=[1,0,0,4;
27    //        5,6,7,8;
28    //        0,0,11,12;
29    //        0,0,0,16];
30    //B=[1,2,0,0]';C=[4,0,0,1];
31    //Sl=ss2ss(syslin('c',A,B,C),rand(A));
32    //[no,X]=contr(Sl(2),Sl(3));CO=X(:,1:no);  //Controllable part
33    //[uo,Y]=unobs(Sl(2),Sl(4));UO=Y(:,1:uo);  //Unobservable part
34    //[Xp,dimc,dimu,dim]=spantwo(CO,UO);    //Kalman decomposition
35    //Slcan=ss2ss(Sl,inv(Xp));
36
37    [X1,dim,dima]=spanplus(A,B);Xp=X1';
38    B1B2B3=Xp*B;B1B2B3=B1B2B3(1:dim,:);
39    [W,dimb]=rowcomp(B1B2B3);W=W(dim:-1:1,:);
40    W11=W(1:dima,1:dima);W21=W(dima+1:dim,1:dima);
41    if rcond(W11)<1.d-10 then
42        //   Which is better?
43        B1B2=B1B2B3(1:dima,:);[W,dimb0]=rowcomp(B1B2);W=W(dima:-1:1,:);
44        [n1,n2]=size(A);
45        Xp=blockdiag(W,eye(n1-dima,n1-dima))*Xp;
46        return;
47    end
48    Q=[eye(dima,dima),zeros(dima,dim-dima);
49    -W21*inv(W11),eye(dim-dima,dim-dima)];
50    Xp(1:dim,:)=Q*W*Xp(1:dim,:);
51endfunction
52