CMSCOPE changed

97 files changed, 2092 insertions, 0 deletions

diff --git a/modules/linear_algebra/macros/aff2ab.bin b/modules/linear_algebra/macros/aff2ab.bin
new file mode 100755
index 000000000..1d0de8e97
--- /dev/null
+++ b/modules/linear_algebra/macros/aff2ab.bin
diff --git a/modules/linear_algebra/macros/aff2ab.sci b/modules/linear_algebra/macros/aff2ab.sci
new file mode 100755
index 000000000..2d86b20e9
--- /dev/null
+++ b/modules/linear_algebra/macros/aff2ab.sci
@@ -0,0 +1,44 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [A,b]=aff2ab(lme,dimX,D,flag)
+    //Y,X,D are lists of matrices.
+    //Y=lme(X,D)= affine fct of Xi's;
+    //[A,b]=matrix representation of lme in canonical basis.
+    // if flag=='sp' A matrix is return in sparse storage.
+    [LHS,RHS]=argn(0)
+    if RHS==3 then flag="f";end
+
+    nvars=0;
+    for k=dimX'
+        nvars=nvars+prod(k);
+    end
+    if part(flag,1)=="f" then
+        x0=zeros(nvars,1);
+        b=list2vec(lme(vec2list(x0,dimX),D));
+        [p,un]=size(b);
+        A=zeros(p,nvars);
+        for k=1:nvars
+            xi=x0;xi(k)=1;
+            A(:,k)=list2vec(lme(vec2list(xi,dimX),D))-b;
+            //   A=[A,sparse(list2vec(lme(vec2list(xi,dimX),D))-b)];
+        end
+    end
+
+    if part(flag,1)=="s" then
+        x0=zeros(nvars,1);
+        b=list2vec(lme(vec2list(x0,dimX),D));
+        A=[];
+        for k=1:nvars
+            xi=x0;xi(k)=1;
+            A=[A,sparse(list2vec(lme(vec2list(xi,dimX),D))-b)];
+        end
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/buildmacros.bat b/modules/linear_algebra/macros/buildmacros.bat
new file mode 100755
index 000000000..a73ca67f8
--- /dev/null
+++ b/modules/linear_algebra/macros/buildmacros.bat
@@ -0,0 +1,11 @@
+
+rem Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+rem Copyright (C) 2008 - INRIA
+rem
+rem This file must be used under the terms of the CeCILL.
+rem This source file is licensed as described in the file COPYING, which
+rem you should have received as part of this distribution.  The terms
+rem are also available at
+rem http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+@..\..\..\bin\scilex -nwni -ns -e exec('buildmacros.sce');quit;
diff --git a/modules/linear_algebra/macros/buildmacros.sce b/modules/linear_algebra/macros/buildmacros.sce
new file mode 100755
index 000000000..9799263d9
--- /dev/null
+++ b/modules/linear_algebra/macros/buildmacros.sce
@@ -0,0 +1,15 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2006-2008 - INRIA - Allan CORNET <allan.cornet@inria.fr>
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+if (isdef("genlib") == %f) then
+    exec(SCI+"/modules/functions/scripts/buildmacros/loadgenlib.sce");
+end
+
+genlib("linear_algebralib","SCI/modules/linear_algebra/macros",%f,%t);
diff --git a/modules/linear_algebra/macros/classmarkov.bin b/modules/linear_algebra/macros/classmarkov.bin
new file mode 100755
index 000000000..4bbf3ec5c
--- /dev/null
+++ b/modules/linear_algebra/macros/classmarkov.bin
diff --git a/modules/linear_algebra/macros/classmarkov.sci b/modules/linear_algebra/macros/classmarkov.sci
new file mode 100755
index 000000000..c68b01b1d
--- /dev/null
+++ b/modules/linear_algebra/macros/classmarkov.sci
@@ -0,0 +1,46 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [perm,rec,tr,indsRec,indsT]=classmarkov(M)
+    //returns a permutation vector perm such that
+    //M(perm,perm) = [M11 0 0 0 0   0]
+    //               [0 M22 0 0     0]
+    //               [0 0 M33       0]
+    //               [      ...      ]
+    //               [0 0       Mrr 0]
+    //               [* *        *  Q]
+    //Each Mii is a Markov matrix of dimension rec(i)  i=1,..,r
+    //Q is sub-Markov matrix of dimension tr
+    //States 1 to sum(rec) are recurrent and states from r+1 to n
+    //are transient.
+    //perm=[indsRec,indsT] where indsRec is a  vector of size sum(rec)
+    //and indsT is a vector of size indsT.
+    if type(M)==1
+        Mb=sparse(M<>0);
+    else Mb=M<>0;
+    end
+    g=mat_2_graph(bool2s(Mb),1,"node-node");
+    [nc,ncomp]=strong_connex(g);
+    indsRec=[];indsT=[];rec=[];tr=0;
+    for i=1:nc
+        inds=find(ncomp==i);
+        nb=size(inds,"*");
+        M1=M(inds,:); M1(:,inds)=[];
+        if sum(M1)==0 then
+            indsRec=[indsRec,inds];
+            rec=[rec,nb];
+        else
+            indsT=[indsT,inds];
+            tr=tr+nb;
+        end
+    end
+    perm=[indsRec,indsT];
+endfunction
+
diff --git a/modules/linear_algebra/macros/cleanmacros.bat b/modules/linear_algebra/macros/cleanmacros.bat
new file mode 100755
index 000000000..4ad1bbae2
--- /dev/null
+++ b/modules/linear_algebra/macros/cleanmacros.bat
@@ -0,0 +1,13 @@
+
+rem Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+rem Copyright (C) 2008 - INRIA
+rem
+rem This file must be used under the terms of the CeCILL.
+rem This source file is licensed as described in the file COPYING, which
+rem you should have received as part of this distribution.  The terms
+rem are also available at
+rem http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+@del *.bin 2>NUL
+@del lib 2>NUL
+@del names 2>NUL
diff --git a/modules/linear_algebra/macros/cmb_lin.bin b/modules/linear_algebra/macros/cmb_lin.bin
new file mode 100755
index 000000000..f6c6b3e69
--- /dev/null
+++ b/modules/linear_algebra/macros/cmb_lin.bin
diff --git a/modules/linear_algebra/macros/cmb_lin.sci b/modules/linear_algebra/macros/cmb_lin.sci
new file mode 100755
index 000000000..812589992
--- /dev/null
+++ b/modules/linear_algebra/macros/cmb_lin.sci
@@ -0,0 +1,19 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [x]=cmb_lin(alfa,x,Beta,y)
+    //  x =cmb_lin(alfa,x,Beta,y) evaluates alfa*x-Beta*y
+    // (low-level routine)
+    //!
+    n=prod(size(x));
+    for j=1:n,
+        x(j)=addf(mulf(alfa,x(j)),mulf(mulf("-1",Beta),y(j)));
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/coff.bin b/modules/linear_algebra/macros/coff.bin
new file mode 100755
index 000000000..9e9ea29ff
--- /dev/null
+++ b/modules/linear_algebra/macros/coff.bin
diff --git a/modules/linear_algebra/macros/coff.sci b/modules/linear_algebra/macros/coff.sci
new file mode 100755
index 000000000..41f6edd01
--- /dev/null
+++ b/modules/linear_algebra/macros/coff.sci
@@ -0,0 +1,33 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [n,d]=coff(m,var)
+    //[N,d]=coff(M [,var]) computes (s*eye-M)^-1
+    // N = numerator polynomial matrix
+    // d = common denominator
+    // var character string ('s' if omitted)
+    // See also : coffg
+    //!
+    //
+    if type(m)<>1 then error(53,1),end
+    if m==[] then n=[];d=1;end
+    [lhs,rhs]=argn(0);if rhs==1 then var="s",end
+    d=clean(poly(m,var)); // denominator
+    [n1,n1]=size(m);
+    for k=1:n1,for l=1:n1,
+            mlk=-m(l,k);
+            if abs(mlk)<1 then mlk=1,end
+            m(l,k)=m(l,k)+mlk;
+            n(k,l)=-(poly(m,var)-d)/mlk;
+            m(l,k)=m(l,k)-mlk
+    end;end
+    if norm(imag(m),1)==0 then n=real(n);d=real(d);end
+    n=clean(n);
+endfunction
diff --git a/modules/linear_algebra/macros/colcomp.bin b/modules/linear_algebra/macros/colcomp.bin
new file mode 100755
index 000000000..16dc4c4b3
--- /dev/null
+++ b/modules/linear_algebra/macros/colcomp.bin
diff --git a/modules/linear_algebra/macros/colcomp.sci b/modules/linear_algebra/macros/colcomp.sci
new file mode 100755
index 000000000..a4dbec36a
--- /dev/null
+++ b/modules/linear_algebra/macros/colcomp.sci
@@ -0,0 +1,36 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [w,rk]=colcomp(a,flag,tol)
+    //Syntaxes : [w,rk]=colcomp(a)
+    //           [w,rk]=colcomp(a,flag)
+    //           [w,rk]=colcomp(a,flag,tol)
+    //
+    //column compression of a i.e. comput. of ker(a)
+    //flag and tol are optional parameters
+    //flag='qr' or 'svd' (defaut 'svd')
+    //tol tolerance parameter (of order %eps as defaut value)
+    //the ma-rk first columns of w span the kernel of a when size(a)=(na,ma)
+
+    [ma,na]=size(a)
+    [lhs,rhs]=argn(0)
+    if a==[] then w=[];rk=0;return;end
+    if norm(a,1) < sqrt(%eps)/10 then rk=0,w=eye(na,na),return,end
+    if rhs ==2 then tol=sqrt(%eps)*norm(a,1)*max(ma,na),end
+    if rhs==1 then flag="svd",tol=sqrt(%eps)*norm(a,1)*max(ma,na);end
+    select flag
+    case "qr" then [q,r,rk,e]=qr(a',tol);
+        //w=[q(:,rk+1:ma),q(:,1:rk)]; <-- le  ma me parait suspect je met na
+        w=q(:,na:-1:1)
+    case "svd" then [u,s,v,rk]=svd(a',tol);
+        //w=[u(:,rk+1:na),u(:,1:rk)];
+        w=u(:,na:-1:1)
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/companion.bin b/modules/linear_algebra/macros/companion.bin
new file mode 100755
index 000000000..20f4026d4
--- /dev/null
+++ b/modules/linear_algebra/macros/companion.bin
diff --git a/modules/linear_algebra/macros/companion.sci b/modules/linear_algebra/macros/companion.sci
new file mode 100755
index 000000000..8934c60da
--- /dev/null
+++ b/modules/linear_algebra/macros/companion.sci
@@ -0,0 +1,50 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function A = companion(p)
+    // Companion matrix.
+    // A=companion(p) is a companion matrix with
+    // characteristic polynomial equal to (or proportional to)
+    // p. If p is a row or column vector of polynomials, the A matrix is block
+    // diagonal and block number i has characteristic polynomial
+    // equal to p(i).
+
+    if type(p) ~= 2
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Polynomial expected."),"companion",1));
+    end
+    // Tranform the row or column vector of poly into a column vector of polynomials
+    p=p(:);
+    // Transpose the vector from column in to row vector,
+    // so that the "for" loop can work properly for each poly.
+    // Caution : ".'", NOT "'"
+    p=p.';
+    A=[];
+    //for pp=p;
+    polynumber = length(p);
+    for polyindex=1:polynumber;
+        pp=p(polyindex);
+        c=coeff(pp);
+        // Reverse the order of the coefficients, so that the coefficient associated with s^n
+        // comes first.
+        c=c($:-1:1);
+        n = length(c);
+        if n <= 1
+            B=[];
+        elseif n == 2
+            B=-c(2)/c(1);
+        else
+            // Caution : ".'", NOT "'"
+            c=c(:).';
+            B=diag(ones(1,n-2),-1);
+            B(1,:) = -c(2:n)/c(1);
+        end
+        A=sysdiag(A,B);
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/cond.bin b/modules/linear_algebra/macros/cond.bin
new file mode 100755
index 000000000..7cc980880
--- /dev/null
+++ b/modules/linear_algebra/macros/cond.bin
diff --git a/modules/linear_algebra/macros/cond.sci b/modules/linear_algebra/macros/cond.sci
new file mode 100755
index 000000000..98f6b8536
--- /dev/null
+++ b/modules/linear_algebra/macros/cond.sci
@@ -0,0 +1,58 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2012- Scilab Enterprises - Adeline CARNIS
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function c = cond(varargin)
+    rhs = argn(2);
+
+    select rhs
+    case 1
+        // c=cond(A) condition number for 2-norm
+        // cond(A) is the ratio of the largest singular of x to the smallest
+        // Warning: the svd function doesn't manage sparse matrices
+        A = varargin(1);
+        if type(A) == 1 then
+            if or(isnan(A)) | or (isinf(A)) then
+                error(msprintf(gettext("%s: Wrong value for input argument #%d: must not contain %s or %s.\n"), "cond", 1, "%nan", "%inf"));
+            end
+            if A == [] then
+                c = 1;
+                return
+            end
+            s = svd(A);
+            if s($) == 0 then
+                c = %inf;
+            else
+                c = s(1)/s($);
+            end
+        else
+            error(msprintf(gettext("%s: Wrong type for input argument #%d: A matrix expected.\n"), "cond", 1));
+        end
+    case 2
+        // c = cond(A, p) where p = 1, %inf, fro, ..
+        // norm(A, p) * norm(inv(A), p)
+        // Warning inv doesn't manage no square and sparse matrices
+        A = varargin(1);
+        p = varargin(2);
+        [m,n] = size(A);
+        if type(A) <> 1 | m <> n then
+            error(msprintf(gettext("%s: Wrong type for input argument #%d: A square matrix expected.\n"),"cond", 1));
+        end
+        if and(type(p) <> [1, 10]) then
+            error(msprintf(gettext("%s: Wrong type for input argument #%d: A scalar or a string expected.\n"),"cond", 2));
+        end
+        if and(p <> [1 2 %inf] & p <> ["inf","fro"]) then
+            error(msprintf(gettext("%s: Wrong value for input argument #%d: must be %d, %d, %s, ''%s'' or ''%s''.\n"), "cond", 2, 1, 2, "%inf", "inf", "fro"));
+        end
+
+        c = norm(A, p) * norm(inv(A), p);
+    else
+        error(msprintf(gettext("%s: Wrong number of input argument(s): %d to %d expected.\n"), "cond", 1, 2));
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/diff.bin b/modules/linear_algebra/macros/diff.bin
new file mode 100755
index 000000000..2c403705b
--- /dev/null
+++ b/modules/linear_algebra/macros/diff.bin
diff --git a/modules/linear_algebra/macros/diff.sci b/modules/linear_algebra/macros/diff.sci
new file mode 100755
index 000000000..a28f1c7ec
--- /dev/null
+++ b/modules/linear_algebra/macros/diff.sci
@@ -0,0 +1,51 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function x=diff(x,N,DIM)
+    rhs=argn(2)
+    if rhs<1 then
+        error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"diff",1,3));
+    end
+    dims=size(x),nd=prod(size(dims))
+    if rhs<3 then DIM="*",end
+    if rhs<2 then N=1,end
+    if DIM=="r" then DIM=1,end
+    if DIM=="c" then DIM=2,end
+    if DIM=="*" then DIM=-1,end
+    if size(DIM,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n"),"diff",1));
+    end
+    DIM=round(DIM)
+    if DIM<>-1&DIM<1 then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Non-negative integers expected.\n"),"diff",3));
+    end
+    if DIM>nd then x=[],return,end
+
+    if type(N)<>1|size(N,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n"),"diff",2));
+    end
+    N=round(N)
+    if N<0 then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Non-negative integers expected.\n"),"diff",3));
+    end
+
+
+    if N==1 then
+        if DIM==-1 then
+            x=x(2:$)-x(1:$-1)
+        else
+            args=emptystr(1,nd)+":"
+            arg1=args;arg1(DIM)="2:$";
+            arg2=args;arg2(DIM)="1:$-1";
+            execstr("x=x("+strcat(arg1,",")+")-x("+strcat(arg2,",")+")")
+        end
+    else
+        for i=1:N, x=diff(x,1,DIM),end
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/eigenmarkov.bin b/modules/linear_algebra/macros/eigenmarkov.bin
new file mode 100755
index 000000000..b07e125c6
--- /dev/null
+++ b/modules/linear_algebra/macros/eigenmarkov.bin
diff --git a/modules/linear_algebra/macros/eigenmarkov.sci b/modules/linear_algebra/macros/eigenmarkov.sci
new file mode 100755
index 000000000..eb0e831ca
--- /dev/null
+++ b/modules/linear_algebra/macros/eigenmarkov.sci
@@ -0,0 +1,49 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [M,Q]=eigenmarkov(P)
+    //Returns normalized left and right eigenvectors
+    //for the eigenvalue 1 of the Markov transition matrix P.
+    //If the multiplicity of this eigenvalue is m and P
+    //is N x N, M is a m x N matrix and Q a N x m matrix.
+    //M(k,:) is the probability distribution vector associated with the kth
+    //ergodic set (recurrent class). M(k,x) is zero if x is not in the
+    //k-th recurrent class.
+    //Q(x,k) is the probability to end in the k-th recurrent class starting
+    //from x.
+    [perm,rec1,tr1,indsRec,indsT]=classmarkov(P)
+    Mn=P(perm,perm);
+    [junk,perminv]=-gsort(-perm);perminv=-perminv;
+    nr=sum(rec1);
+    T=Mn(nr+1:$,nr+1:$);L=Mn(nr+1:$,1:nr);
+    p=0;V=[];
+    for k=rec1
+        v=L(:,p+1:p+k);V=[V,sum(v,"c")];
+        p=p+k;
+    end
+    LL=zeros(nr,size(rec1,"*"));
+    p=0;
+    for k=1:size(rec1,"*")
+        LL(p+1:p+rec1(k),k)=1;
+        p=p+rec1(k);
+    end
+    Q=[LL;inv(eye()-T)*V];
+    Q=Q(perminv,:);
+    p=0;M=[];
+    for kk=1:size(rec1,"*")
+        classe=p+1:p+rec1(kk);
+        p=p+rec1(kk);
+        Mres=Mn(classe,classe);
+        w=kernel((Mres-eye())')';
+        M=sysdiag(M,w./sum(w));
+    end
+    M=[M,zeros(size(M,1),size(P,1)-size(M,2))];
+    M=M(:,perminv);
+endfunction
diff --git a/modules/linear_algebra/macros/fullrf.bin b/modules/linear_algebra/macros/fullrf.bin
new file mode 100755
index 000000000..fa79933b9
--- /dev/null
+++ b/modules/linear_algebra/macros/fullrf.bin
diff --git a/modules/linear_algebra/macros/fullrf.sci b/modules/linear_algebra/macros/fullrf.sci
new file mode 100755
index 000000000..023461142
--- /dev/null
+++ b/modules/linear_algebra/macros/fullrf.sci
@@ -0,0 +1,30 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Q,M,rk]=fullrf(A,tol)
+    //[Q,M,rk]=fullrf(A)
+    //Full rank factorization : A=Q.M
+    //with range(Q)=range(A) and ker(M)=ker(A),
+    //Q full column rank , M full row rank
+    // rk = rank(A) = #columns(Q) = #rows(M)
+    //F.D.
+    //!
+    [lhs,rhs]=argn(0)
+    na1=norm(A,1);
+    if rhs==1 then tol=sqrt(%eps);end
+    if na1 < 1.d-10 then Q=[];M=[];rk=0;return;end
+    tol1=tol*na1;
+    [U,s,V,rk]=svd(A,tol1);
+    sq=sqrt(s);
+    Q=U*sq;M=sq*V';
+    Q=Q(:,1:rk);M=M(1:rk,:);
+endfunction
+
+
diff --git a/modules/linear_algebra/macros/fullrfk.bin b/modules/linear_algebra/macros/fullrfk.bin
new file mode 100755
index 000000000..2290c78a8
--- /dev/null
+++ b/modules/linear_algebra/macros/fullrfk.bin
diff --git a/modules/linear_algebra/macros/fullrfk.sci b/modules/linear_algebra/macros/fullrfk.sci
new file mode 100755
index 000000000..4ae71ec9d
--- /dev/null
+++ b/modules/linear_algebra/macros/fullrfk.sci
@@ -0,0 +1,38 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Bk,Ck]=fullrfk(A,k)
+    // This macro computes the full rank factorization of A^k i.e.
+    // Bk*Ck=A^k where Bk is full column rank and Ck full row rank.
+    // One has Range(Bk)=Range(A^k) and ker(Ck)=ker(A^k).
+    // For k=1 fullrfk is the same as fullrf.
+    // F.D (1990)
+    //!
+    [lhs,rhs]=argn(0);
+    if rhs==1,
+        k=1;
+    end
+    [n,n]=size(A);
+    //k=0
+    if k==0,
+        Bk=eye(n,n);Ck=Bk;
+        return;
+    end
+    //k=1
+    if k==1,
+        [Bk,Ck]=fullrf(A);
+        return;
+    end
+    [Bk,Ck]=fullrf(A);B=Bk;C=Ck;
+    for l=2:k
+        [B,C,dim]=fullrf(C*B);
+        Bk=Bk*B;Ck=C*Ck;     // Bk*Ck = A^k  (Full rank factorization)
+    end;
+endfunction
diff --git a/modules/linear_algebra/macros/genmarkov.bin b/modules/linear_algebra/macros/genmarkov.bin
new file mode 100755
index 000000000..65072b774
--- /dev/null
+++ b/modules/linear_algebra/macros/genmarkov.bin
diff --git a/modules/linear_algebra/macros/genmarkov.sci b/modules/linear_algebra/macros/genmarkov.sci
new file mode 100755
index 000000000..cecb7797c
--- /dev/null
+++ b/modules/linear_algebra/macros/genmarkov.sci
@@ -0,0 +1,62 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function M=genmarkov(rec,tr,flag)
+    //Returns in M a random Markov transition probability matrix
+    //with size(rec,1) recurrent classes with rec(1),...rec($)
+    //entries respectively and tr transient states.
+    // If the optional parameter flag='perm' is entered a random
+    //permutation of the states is performed.
+    [lhs,rhs]=argn(0);
+    if rhs==2 then flag="noperm";end
+    M=[];r=sum(rec);
+    for k=rec
+        m=rand(k,k,"u");
+        m=m./(sum(m,"c")*ones(1,k));
+        M=sysdiag(M,m);
+    end
+    if type(tr) ~= 15 then
+        n=r+tr;
+        MT=rand(tr,n,"u");MT=MT./(sum(MT,"c")*ones(1,n));
+        M=[[M,zeros(r,tr)];MT];
+        if flag=="perm" then [p,q]=gsort(rand(n,1,"u"));M=M(q,q);end
+    else
+        // tr=list(n1,[a1,a2,...ar],n2,[b1,...br],...)
+        l=size(tr)/2;   //2*size(rec,2)
+        Q=[];
+        for kk=1:l
+            nt=tr(1+2*(kk-1));
+            Q=sysdiag(Q,rand(nt,nt));
+        end
+        Nt=size(Q,1);
+        L=[];
+        nclrec=size(rec,"*");
+        for kk=1:l
+            L1=[];indi=tr(2+2*(kk-1));nt=tr(1+2*(kk-1));
+            for i=1:nclrec
+                if indi(i)==0 then L1=[L1,zeros(nt,rec(i))]; else L1=[L1,rand(nt,rec(i))];end
+            end
+            L=[L;L1]
+        end
+        LQ=[L,Q];
+        LQ=LQ./(sum(LQ,"c")*ones(1,size(LQ,2)));
+        M=[[M,zeros(size(M,1),size(Q,2))];LQ];
+    end
+    if flag=="perm" then [p,q]=gsort(rand(size(M,1),1,"u"));M=M(q,q);end
+
+    //n=size(M,1);[p,q]=gsort(rand(n,1));M1=M(q,q);M=M1
+
+    //Ms=sparse(M);Ms(find(Ms~=0))=1;
+    //G=mat_2_graph(Ms,1,'node-node');show_graph(G);
+    //[nc,ncomp]=strong_connex(G)
+    //[pp,qq]=gsort(ncomp);Mg=M(qq,qq)
+    //Mgs=sparse(Mg);
+    //Mgs(find(Mgs~=0))=1;Gg=mat_2_graph(Mgs,1,'node-node');show_graph(Gg);
+endfunction
diff --git a/modules/linear_algebra/macros/givens.bin b/modules/linear_algebra/macros/givens.bin
new file mode 100755
index 000000000..cd4a6ed93
--- /dev/null
+++ b/modules/linear_algebra/macros/givens.bin
diff --git a/modules/linear_algebra/macros/givens.sci b/modules/linear_algebra/macros/givens.sci
new file mode 100755
index 000000000..8c0441f15
--- /dev/null
+++ b/modules/linear_algebra/macros/givens.sci
@@ -0,0 +1,34 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [u,c]=givens(x,y)
+    //Syntax : u=givens(xy)
+    //         u=givens(x,y)
+    //
+    // xy = [x;y], u=givens(xy)
+    // returns a 2*2 matrix u such that u*xy=[r;0].
+    // c is equal to u*xy
+    // givens(x,y)=givens([x;y])
+    //
+    //!
+    [lhs,rhs]=argn(0);
+    if rhs==2 then x=[x;y];end
+    if or(size(x)<>[2 1]) then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A column vector expected.\n"),"givens",1));
+    end
+    if x(2)<>0 then
+        r = norm(x);
+        u = [x'; -x(2) x(1)]/r;
+        c = [r; 0];
+    else
+        u=eye(2,2)
+        c=x
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/glever.bin b/modules/linear_algebra/macros/glever.bin
new file mode 100755
index 000000000..8c8299d0d
--- /dev/null
+++ b/modules/linear_algebra/macros/glever.bin
diff --git a/modules/linear_algebra/macros/glever.sci b/modules/linear_algebra/macros/glever.sci
new file mode 100755
index 000000000..135d96f31
--- /dev/null
+++ b/modules/linear_algebra/macros/glever.sci
@@ -0,0 +1,52 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 1988-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Bfs,Bis,chis]=glever(E,A,s)
+    // [Bfs,Bis,chis]=glever(E,A,'s')
+    // Computation of (s*E-A)^-1 ('s'=character string with default value 's')
+    // Generelized Leverrier's algorithm for a matrix pencil ;
+    // (s*E-A)^-1 = (Bfs/chis) - Bis
+    // chis = characteristic polynomial (up to a multiplicative constant)
+    // Bfs  = polynomial matrix
+    // Bis  = polynomial matrix ( - expansion of (s*E-A)^-1 at infinity).
+    // Caveat: uses clean to simplify Bfs,Bis and chis !
+    // See also shuffle, determ, invr, coffg
+
+    [LHS,RHS]=argn(0);
+    if RHS==1 then [E,A]=pen2ea(E),s=poly(0,"s");end
+    if RHS==2 then s=poly(0,"s"),end;
+    if RHS==3 then s=poly(0,s);end
+    [Si,Pi,Di,index]=penlaur(E,A);
+    k=round(sum(diag(Si*E)));
+
+    a0=1;
+    B=Si;
+    SiASi=Si*A*Si;
+    chis=a0+0*s;
+    Bfs=Si+0*s*Si;
+
+    for i=1:k,
+        B=SiASi*E*B;
+        alfa=-sum(diag(E*B))/i;
+        B=B+alfa*Si;
+        chis=s*chis+alfa;
+        if i<k then Bfs=s*Bfs+B;end
+    end
+    Bis=Pi;
+    AAD=s*A*Di;
+    P=eye(A);
+
+    for nu=1:index+1,
+        P=AAD*P;
+        Bis=clean(Bis+Pi*P,1.d-10);
+    end
+    Bfs=clean(Bfs,1.d-10);
+    chis=clean(chis,1.d-10);
+endfunction
diff --git a/modules/linear_algebra/macros/gschur.bin b/modules/linear_algebra/macros/gschur.bin
new file mode 100755
index 000000000..d3bfa8286
--- /dev/null
+++ b/modules/linear_algebra/macros/gschur.bin
diff --git a/modules/linear_algebra/macros/gschur.sci b/modules/linear_algebra/macros/gschur.sci
new file mode 100755
index 000000000..25c23b79e
--- /dev/null
+++ b/modules/linear_algebra/macros/gschur.sci
@@ -0,0 +1,47 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [As,Es,Z,dim]= gschur(A,E,extern)
+
+    if argn(2)<=2 then
+        warning("Obsolete function. Please use schur instead.")
+        if argn(1)==2 then
+            [As,Es]= schur(A,E)
+        else
+            [As,Es,Q,Z]= schur(A,E);Q=Q'
+        end
+    elseif type(extern)==10 then
+        if convstr(extern)=="c"|convstr(extern)=="d" then
+            warning(msprintf(gettext("%s: Obsolete function. Please replace ''%s'' by ''%s''."),"gschur","gschur","schur"));
+            if argn(1)==4 then
+                [As,Es,Z,dim]= schur(A,E,extern)
+            elseif argn(1)==2 then
+                [As,Es,Z,dim]= schur(A,E,extern)
+                Es=dim;As=Z;
+            end
+        else // hard coded
+            warning(msprintf(gettext("%s: Obsolete function. Former external functions cannot be used anymore."),"gschur"));
+            //impossible to redefine
+        end
+    else //coded by a scilab function
+        //---- old------------------
+        //flag=extern(x)
+        //x(1) ==1 ==> x(2:3)=[al,be]
+        //x(1) ==2 ==> x(2:3)=[s,p]
+        warning(msprintf(gettext("%s: Obsolete function. Please replace ''%s'' by ''%s''."),"gschur","gschur","schur"));
+        deff("t=%_rule(Alpha,Beta)",["if imag(Alpha)==0 then"
+        "  t=extern([1,real(Alpha),Beta])==1"
+        "else"
+        "  c=Alpha/Beta"
+        "  t=extern([2,real(c+c''),real(c*c'')])==1"
+        "end;"])
+        [As,Es,Z,dim]= schur(A,E,%_rule)
+    end
+endfunction
+
diff --git a/modules/linear_algebra/macros/gspec.bin b/modules/linear_algebra/macros/gspec.bin
new file mode 100755
index 000000000..743b453f9
--- /dev/null
+++ b/modules/linear_algebra/macros/gspec.bin
diff --git a/modules/linear_algebra/macros/gspec.sci b/modules/linear_algebra/macros/gspec.sci
new file mode 100755
index 000000000..275109ef4
--- /dev/null
+++ b/modules/linear_algebra/macros/gspec.sci
@@ -0,0 +1,15 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [al,be,Z]=gspec(A,B)
+    //For backward compatibility
+    warning(msprintf(gettext("%s: Obsolete function. Please replace ''%s'' by ''%s''."),"gspec","gspec","spec"));
+    [al,be,Z]=spec(A,B)
+endfunction
diff --git a/modules/linear_algebra/macros/householder.bin b/modules/linear_algebra/macros/householder.bin
new file mode 100755
index 000000000..1ffbdbb80
--- /dev/null
+++ b/modules/linear_algebra/macros/householder.bin
diff --git a/modules/linear_algebra/macros/householder.sci b/modules/linear_algebra/macros/householder.sci
new file mode 100755
index 000000000..f25e05832
--- /dev/null
+++ b/modules/linear_algebra/macros/householder.sci
@@ -0,0 +1,25 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function u=householder(v,w)
+    //Syntax
+    //u=householder(v [,w])
+    //Description
+    //given 2 column vectors v w of same size householder(v,w) returns a unitary
+    //column vector u, such that (eye-2*u*u')*v is proportional to w.
+    //(eye-2*u*u') is the orthogonal Householder reflexion matrix
+    //
+    // w default value is eye(v). In this case vector (eye-2*u*u')*v is the
+    // vector  eye(v)*(+-norm(v))
+    [lhs,rhs]=argn(0)
+    if rhs<2 then w=eye(v),end
+    a=-sqrt((v'*v)/(w'*w))
+    u=v+a*w
+    u=u/norm(u)
+endfunction
diff --git a/modules/linear_algebra/macros/im_inv.bin b/modules/linear_algebra/macros/im_inv.bin
new file mode 100755
index 000000000..56924abac
--- /dev/null
+++ b/modules/linear_algebra/macros/im_inv.bin
diff --git a/modules/linear_algebra/macros/im_inv.sci b/modules/linear_algebra/macros/im_inv.sci
new file mode 100755
index 000000000..99ea513d2
--- /dev/null
+++ b/modules/linear_algebra/macros/im_inv.sci
@@ -0,0 +1,45 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [X,dim,Y]=im_inv(A,B,tol)
+    //[X,dim]=im_inv(A,B [,tol]) computes (A^-1)(B) i.e vectors whose
+    // image through A are in range(B).
+    // The dim first columns de X span (A^-1) (B)
+    // tol is a threshold to test if a  subspace is included in an other
+    // default value tol = 100*%eps;
+
+    [lhs,rhs]=argn(0);
+    [nA,mA]=size(A);[nB,mB]=size(B);
+    if rhs==2 then tol=100*%eps*mA*nA*nB*mB,end;
+    if nA<>nB then
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"),"im_inv",1,2))
+        return
+    end
+    // basis for im(B)
+    [Y,rB]=rowcomp(B);//u=Y'
+
+    //Trivial cases
+    if rB >= nA then X=eye(mA,mA);dim=mA,return;end;
+    if rB ==0 then [X,k1]=colcomp(A);dim=mA-k1,return,end
+    //
+    up=1:rB;low=rB+1:nA;
+    A=Y*A;   //update
+
+    //vectors with image in B
+    [X,r1]=colcomp(A(low,:))
+    A1=A*X;    //update
+    Aup=A1(up,:);
+    Alow=A1(low,:);    //Alow(:,1:ma-r1)=0 by construction
+    if norm(Alow,1) <= tol*norm(Aup,1) then dim=mA,return,end
+    dim=mA-r1;
+endfunction
+
+
+
diff --git a/modules/linear_algebra/macros/kernel.bin b/modules/linear_algebra/macros/kernel.bin
new file mode 100755
index 000000000..9b4a7110f
--- /dev/null
+++ b/modules/linear_algebra/macros/kernel.bin
diff --git a/modules/linear_algebra/macros/kernel.sci b/modules/linear_algebra/macros/kernel.sci
new file mode 100755
index 000000000..931492acc
--- /dev/null
+++ b/modules/linear_algebra/macros/kernel.sci
@@ -0,0 +1,65 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function kerA=kernel(A,tol,flag)
+    //Syntax : [x]=kernel(A [,tol [,flag]])
+    // flag = 'svd' or 'sqr' no flag if A sparse
+    if type(A)==1 then
+        [lhs,rhs]=argn(0);
+        [n,m]=size(A);
+        select rhs
+        case 1 then
+            [kerA,rk]=colcomp(A);
+        case 2 then
+            [kerA,rk]=colcomp(A,"svd",tol)
+        case 3 then
+            [kerA,rk]=colcomp(A,flag,tol)
+        end
+        kerA=kerA(:,1:m-rk);return;
+    end
+    if type(A)==5 then
+        [lhs,rhs]=argn(0);
+        [ma,na]=size(A);
+        if rhs==3 then
+            error(msprintf(gettext("%s: SVD and QR not implemented in sparse.\n"),"kernel"));
+        end
+        if rhs==2 then tol=1.d-10;else tol=%eps;end
+        %tol=tol*max(abs(A))*max(ma,na);
+        if ma<na then
+            A=[A;sparse([],[],[na-ma,na])];
+        end
+        if ma>na then
+            A=A'*A;ma=na;
+        end
+        [ptr,rkA]=lufact(A,[%tol,0.001]);
+        [P,L,U,Q]=luget(ptr);
+        ludel(ptr);
+        U=clean(U);
+        [ptrU,rk]=lufact(U,[%tol,.001]);
+        nma=max(na,ma);
+        Y=[sparse([],[],[rkA,nma-rkA]);speye(nma-rkA,nma-rkA)];
+        kerA=[];
+        for k=1:na-rkA
+            bb=full(Y(:,k));
+            ww=sparse(lusolve(ptrU,bb));
+            kerA=[kerA,ww];
+        end
+        ludel(ptrU);
+        if na<>rkA then
+            kerA=clean(Q'*kerA);
+        end
+        if ma>na then kerA=kerA(1:na,:);
+        end
+    else
+        error(msprintf(gettext("%s: This feature has not been implemented.\n"),"kernel"));
+    end
+endfunction
+
+
diff --git a/modules/linear_algebra/macros/kroneck.bin b/modules/linear_algebra/macros/kroneck.bin
new file mode 100755
index 000000000..d2dbe9440
--- /dev/null
+++ b/modules/linear_algebra/macros/kroneck.bin
diff --git a/modules/linear_algebra/macros/kroneck.sci b/modules/linear_algebra/macros/kroneck.sci
new file mode 100755
index 000000000..472afb479
--- /dev/null
+++ b/modules/linear_algebra/macros/kroneck.sci
@@ -0,0 +1,67 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Q,Z,Qd,Zd,numbeps,numbeta]=kroneck(E,A)
+    //  Kronecker form:
+    //
+    //            | sE(eps)-A(eps) |        X       |      X     |      X        |
+    //            |----------------|----------------|------------|---------------|
+    //            |        O       | sE(inf)-A(inf) |      X     |      X        |
+    // Q(sE-A)Z = |---------------------------------|----------------------------|
+    //            |                |                |            |               |
+    //            |        0       |       0        | sE(f)-A(f) |      X        |
+    //            |--------------------------------------------------------------|
+    //            |                |                |            |               |
+    //            |        0       |       0        |      0     | sE(eta)-A(eta)|
+    //
+    //  Ec=Q*E*Z, Ac=Q*A*Z, eps=Qd(1) x Zd(1) ,inf=Qd(2) x Zd(2)
+    //  f = Qd(3) x Zd(3),  eta=Qd(4)xZd(4)
+    //
+    // numbeps(1) = # of eps blocks of size 0 x 1
+    // numbeps(2) = # of eps blocks of size 1 x 2
+    // numbeps(3) = # of eps blocks of size 2 x 3     etc...
+    //
+    // numbeta(1) = # of eta blocks of size 1 x 0
+    // numbeta(2) = # of eta blocks of size 2 x 1
+    // numbeta(3) = # of eta blocks of size 3 x 2     etc...
+    //
+    // interface  F.D. from Slicot-fstair
+    // T. Beelen's routines
+    [LHS,RHS]=argn(0);
+    if RHS==1 then [E,A]=pen2ea(E);end
+    [Q,Z,Ec,Ac,Qd,Zd,numbeps]=quaskro(E,A);
+    rows=Qd(1)+Qd(2)+1:Qd(1)+Qd(2)+Qd(3);
+    cols=Zd(1)+Zd(2)+1:Zd(1)+Zd(2)+Zd(3);
+    if cols==[] then
+        Zd(4)=Zd(3);Qd(4)=Qd(3);Qd(3)=0;
+        numbeta=Qd(4)-Zd(4);
+        return;
+    end;
+    if rows==[] then Qd(4)=0;Zd(4)=0;numbeta=0;return;end;
+    Er=Ec(rows,cols);
+    Ar=Ac(rows,cols);
+    E1=pertrans(Er);
+    A1=pertrans(Ar);
+    [Q2,Z2,Ec2,Ac2,Qd2,Zd2,numbeta]=quaskro(E1,A1);
+    // Here check that...
+    //Zd2(2)=0  Qd2(2)=0;
+    //Zd2(1)+Zd2(3)=Zd(3)
+    //Qd2(1)+Qd2(3)=Qd(3)
+    Z1=pertrans(Q2);
+    Q1=pertrans(Z2);
+    Zd1=[Qd2(3),Qd2(1)];
+    Qd1=[Zd2(3),Zd2(1)];
+    Zd=[Zd(1),Zd(2),Zd1];
+    Qd=[Qd(1),Qd(2),Qd1];
+    Z(:,cols)=Z(:,cols)*Z1;
+    Q(rows,:)=Q1*Q(rows,:);
+endfunction
+
+
diff --git a/modules/linear_algebra/macros/lib b/modules/linear_algebra/macros/lib
new file mode 100755
index 000000000..7160a333c
--- /dev/null
+++ b/modules/linear_algebra/macros/lib
diff --git a/modules/linear_algebra/macros/linsolve.bin b/modules/linear_algebra/macros/linsolve.bin
new file mode 100755
index 000000000..031d88f81
--- /dev/null
+++ b/modules/linear_algebra/macros/linsolve.bin
diff --git a/modules/linear_algebra/macros/linsolve.sci b/modules/linear_algebra/macros/linsolve.sci
new file mode 100755
index 000000000..1e9a11cf3
--- /dev/null
+++ b/modules/linear_algebra/macros/linsolve.sci
@@ -0,0 +1,77 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [x0,kerA]=linsolve(A,b,x0)
+    %tol=1.D-10;
+    // Finds all x solution to Ax+b=0;
+    // x0=particular solution; kerA=nullspace of A
+    // Any x=x0+kerA*w with arbitrary w solves A*x+b=0
+    [LHS,RHS]=argn(0);
+    select type(A)
+    case 1 then        //full matrix
+        Ab=[A,b];[ma,na]=size(Ab);
+        [W,rk]=colcomp(Ab);
+        W=W(:,1:na-rk);last=W(na,:);
+        [W2,rk1]=colcomp(last);
+        if rk1==0 then
+            warning(gettext("Conflicting linear constraints."));
+            x0=[];kerA=[];return;
+        end
+        W=W*W2;
+        kerA=W(1:na-1,1:na-rk-1);
+        if RHS==3 then
+            if norm(A*x0+b,1)<%tol then
+                return;
+            end
+            disp(gettext("Recomputing initial guess"));
+        end
+        piv=W(na,na-rk);x0=W(1:na-1,na-rk)/piv;
+    case 5 then        //sparse matrix
+
+        [ma,na]=size(A);
+        %tol=1.D-10*max(abs(A))*max(ma,na);
+        if ma<na then
+        A=[A;sparse([],[],[na-ma,na])];b=[b;zeros(na-ma,1)];end
+        if ma>na then
+            //A=[A,sparse([],[],[ma,ma-na])];x0=[x0;zeros(ma-na,1)];
+            b=A'*b;A=A'*A;ma=na;
+        end
+        [ptr,rkA]=lufact(A,[%tol,0.001]);
+        [P,L,U,Q]=luget(ptr);
+        if RHS==3 then
+            if norm(A*x0+b,1)>%tol then
+                x0=lusolve(ptr,-b);
+                disp("recomputing initial guess");
+            end
+        end
+        if RHS==2 then
+            x0=lusolve(ptr,-b);
+        end
+        ludel(ptr);
+        U=clean(U);
+        [ptrU,rk]=lufact(U,[%tol,.001]);
+        nma=max(na,ma);
+        Y=[sparse([],[],[rkA,nma-rkA]);speye(nma-rkA,nma-rkA)];
+        kerA=[];
+        for k=1:na-rkA
+            bb=full(Y(:,k));
+            ww=sparse(lusolve(ptrU,bb));
+            kerA=[kerA,ww];
+        end
+        if na<>rkA then
+            kerA=clean(Q'*kerA);
+        end
+        if norm(A*x0+b,1)>%tol then
+            warning(msprintf(gettext("Possible Conflicting linear constraints, error in the order of %s"),string(norm(A*x0+b,1)) ));
+        end
+        if ma>na then kerA=kerA(1:na,:);x0=x0(1:na,1);end
+        ludel(ptrU);
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/names b/modules/linear_algebra/macros/names
new file mode 100755
index 000000000..fdee597e3
--- /dev/null
+++ b/modules/linear_algebra/macros/names
@@ -0,0 +1,46 @@
+aff2ab
+classmarkov
+cmb_lin
+coff
+colcomp
+companion
+cond
+diff
+eigenmarkov
+fullrf
+fullrfk
+genmarkov
+givens
+glever
+gschur
+gspec
+householder
+im_inv
+kernel
+kroneck
+linsolve
+nlev
+orth
+pbig
+pen2ea
+pencan
+penlaur
+pinv
+polar
+proj
+projspec
+psmall
+quaskro
+randpencil
+range
+rank
+rowcomp
+rowshuff
+rref
+spaninter
+spanplus
+spantwo
+sqroot
+squeeze
+sva
+trace
diff --git a/modules/linear_algebra/macros/nlev.bin b/modules/linear_algebra/macros/nlev.bin
new file mode 100755
index 000000000..5c6824cac
--- /dev/null
+++ b/modules/linear_algebra/macros/nlev.bin
diff --git a/modules/linear_algebra/macros/nlev.sci b/modules/linear_algebra/macros/nlev.sci
new file mode 100755
index 000000000..9a2fe3921
--- /dev/null
+++ b/modules/linear_algebra/macros/nlev.sci
@@ -0,0 +1,65 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 1987-2008 - INRIA - Serge STEER <serge.steer@inria.fr>
+// Copyright (C) 1987-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function[m,den]=nlev(a,z,rmax)
+    //[num,den]=nlev(a,z [,rmax])  calcule (z*eye-a)**(-1) par une version
+    //modifiee de l'algorithme de leverrier numeriquement plus stable.
+    //
+    //z     chaine de caracteres donnant le nom de la variable formelle
+    //rmax  parametre optionnel utilise pour bloc diagonaliser a (voir la
+    //      fonction bdiag)
+    //
+
+    rhs=argn(2);
+    if rhs < 2 then
+        error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"),"nlev",2,3));
+    end
+
+    z=poly(0,z);
+
+    if rhs==3 then
+        [a,x,bs]=bdiag(a,rmax),
+    else
+        [a,x,bs]=bdiag(a),
+    end
+
+    [m1 n1]=size(a)
+    if m1<>n1 then
+        error(20,1);
+    end
+    k=1;
+    m=[];
+    v=ones(1,n1);
+    den=1;
+    for n=bs';k1=k:k-1+n;
+        //       algorithme de leverrier
+        h=z*eye(n,n)-a(k1,k1)
+        f=eye(n,n)
+        for kl=1:n-1,
+            b=h*f,
+            d=-sum(diag(b))/kl,
+            f=b+eye()*d,
+        end
+        d=sum(diag(h*f))/n
+        //
+        den=den*d;
+        l=[1:k-1,k+n:n1] ,
+        if l<>[] then
+            v(l)=v(l)*d;
+        end
+        m=[m,x(:,k1)*f];
+        k=k+n;
+    end;
+    m=m*diag(v)*inv(x);
+endfunction
+
+
+
diff --git a/modules/linear_algebra/macros/orth.bin b/modules/linear_algebra/macros/orth.bin
new file mode 100755
index 000000000..eb58dab1f
--- /dev/null
+++ b/modules/linear_algebra/macros/orth.bin
diff --git a/modules/linear_algebra/macros/orth.sci b/modules/linear_algebra/macros/orth.sci
new file mode 100755
index 000000000..9fd49fd40
--- /dev/null
+++ b/modules/linear_algebra/macros/orth.sci
@@ -0,0 +1,18 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Q]=orth(A)
+    // Orthogonal basis for the span of A.
+    // Range(Q) = Range(A) and Q'*Q=eye
+    //!
+    [X,n]=rowcomp(A);X=X';
+    Q=X(:,1:n);
+endfunction
+
diff --git a/modules/linear_algebra/macros/pbig.bin b/modules/linear_algebra/macros/pbig.bin
new file mode 100755
index 000000000..9ae7d96c3
--- /dev/null
+++ b/modules/linear_algebra/macros/pbig.bin
diff --git a/modules/linear_algebra/macros/pbig.sci b/modules/linear_algebra/macros/pbig.sci
new file mode 100755
index 000000000..a9e78fca3
--- /dev/null
+++ b/modules/linear_algebra/macros/pbig.sci
@@ -0,0 +1,49 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Q,M]=pbig(A,thres,flag)
+    //Projection on eigensubspace associated with eigenvalues
+    //with real part >= thres (flag='c') or with modulus >= thres (flag='d')
+    //Projection is defined by Q*M. Eigenvalues of M*A*Q = eigenvalues
+    //of A with real part >= thres (case 'c',...).
+    //If [Q1,M1]== full rank factorization (fullrf) of eye-Q*M then evals of
+    // M1*A*Q1 = evals of A with real part < thres (case 'c',...).
+    // See also psmall.
+    //!
+
+    [n,n]=size(A);
+    thres=real(thres);
+    if flag=="c" then
+        deff("[flag]=%smallei(x)","flag=real(x) >= thres")
+        deff("[flag]=%bigeig(x)","flag=real(x) < thres")
+    elseif flag=="d" then
+        deff("[flag]=%smallei(x)","flag=abs(x) >= thres")
+        deff("[flag]=%bigeig(x)","flag=abs(x) < thres")
+    else
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"pbig",3,"''c'',''d''"));
+    end
+    //
+    [X,dsmall] = schur(A,%smallei);
+    [Y,dbig] = schur(A,%bigeig);
+    Q=X(:,1:dsmall);
+    if Q==[] then M=[];return;end
+    Y1=Y';
+    M1=Y1(dbig+1:n,:);
+    E=M1*Q;
+    if rcond(E)>1.d-7 then
+        M=E\M1;
+    else
+        [Ab,X0]=balanc(A);
+        [X,dsmall] = schur(Ab,%smallei);X1=X*X0;Q=X1(:,1:dsmall);
+        [Y,dbig] = schur(Ab,%bigeig);Y1=inv(X0)*Y';M=Y1(dbig+1:n,:);
+        E=M*Q;
+        M=E\M;
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/pen2ea.bin b/modules/linear_algebra/macros/pen2ea.bin
new file mode 100755
index 000000000..a4bf4499e
--- /dev/null
+++ b/modules/linear_algebra/macros/pen2ea.bin
diff --git a/modules/linear_algebra/macros/pen2ea.sci b/modules/linear_algebra/macros/pen2ea.sci
new file mode 100755
index 000000000..17dbf0ede
--- /dev/null
+++ b/modules/linear_algebra/macros/pen2ea.sci
@@ -0,0 +1,14 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [e,a]=pen2ea(f)
+    e=coeff(f,1);
+    a=-coeff(f,0);
+endfunction
diff --git a/modules/linear_algebra/macros/pencan.bin b/modules/linear_algebra/macros/pencan.bin
new file mode 100755
index 000000000..0431f3dc5
--- /dev/null
+++ b/modules/linear_algebra/macros/pencan.bin
diff --git a/modules/linear_algebra/macros/pencan.sci b/modules/linear_algebra/macros/pencan.sci
new file mode 100755
index 000000000..77b551000
--- /dev/null
+++ b/modules/linear_algebra/macros/pencan.sci
@@ -0,0 +1,30 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Q,M,i1]=pencan(E,A)
+    // [Q,M,i1]=pencan(E,A)
+    // Given the pencil s*E-A pencan returns matrices Q and M
+    // such than M*(s*E-A)*Q is in "canonical" form i.e.
+    // M*E*Q is a block matrix [I,0;
+    //                         0,N]    with N nilpotent
+    // and i1 = size of the I matrix above (# of finite zeros of (sE-A)).
+    // M*A*Q is a block matrix [Ar,0;
+    //                         0,I ]
+    //See glever,  penlaur
+    //!
+    [LHS,RHS]=argn(0);
+    if RHS==1 then [E,A]=pen2ea(E);end
+    [Si,Pi,Di,index]=penlaur(E,A);
+    [Q1,M1]=fullrf(Si);
+    [Q2,M2]=fullrf(Pi);
+    [i1,i2]=size(M1);
+    M=[M1;M2];
+    Q=[Q1,Q2];
+endfunction
diff --git a/modules/linear_algebra/macros/penlaur.bin b/modules/linear_algebra/macros/penlaur.bin
new file mode 100755
index 000000000..9d24cf952
--- /dev/null
+++ b/modules/linear_algebra/macros/penlaur.bin
diff --git a/modules/linear_algebra/macros/penlaur.sci b/modules/linear_algebra/macros/penlaur.sci
new file mode 100755
index 000000000..27d57d372
--- /dev/null
+++ b/modules/linear_algebra/macros/penlaur.sci
@@ -0,0 +1,43 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 1988-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Si,Pi,Di,order]=penlaur(E,A)
+    //[Si,Pi,Di,order]=penlaur(E,A)
+    // First Laurent coefficients of (s*E-A)^-1;
+    // (s*E-A)^-1 = ... + Si/s - (Pi + s*Di + ... +s^order Ni) at s = infinity
+    // order = order of the singularity
+    // The matrix s*E-A should be invertible.
+    // Experimental version: troubles when bad conditioning of
+    // (so*E-A)...)
+    //!
+    [LHS,RHS]=argn(0);
+    if RHS==1 then [E,A]=pen2ea(E);end
+    seed=rand("seed");typ=rand("info");
+    rand("normal");rand("seed",0);
+    tests=rand(1,10);
+    conditions=0*tests;k=1;
+    for s0=tests, conditions(k)=cond(s0*E-A);k=k+1;end
+    [w,k1]=min(conditions);
+    rand(typ);rand("seed",seed)
+    if w>1.d+20 then
+        error(msprintf(gettext("%s: Singular pencil."),"penlaur"));
+        return;
+    end
+    s0=tests(k1);
+    J=inv(s0*E-A);
+    [Se,Pe,De,i1]=projspec(J*E);
+    [Sa,Pa,Da,i2]=projspec(J*A);
+    order=i1-1;
+    Si=Se*J;
+    Pi=Pe*Sa*J;
+    Di=Pi*E*Pi;
+    if order==0 then Di=0*Di;end
+    //[S,P,D,index]=projspec(Di*E);
+endfunction
diff --git a/modules/linear_algebra/macros/pinv.bin b/modules/linear_algebra/macros/pinv.bin
new file mode 100755
index 000000000..81b456f14
--- /dev/null
+++ b/modules/linear_algebra/macros/pinv.bin
diff --git a/modules/linear_algebra/macros/pinv.sci b/modules/linear_algebra/macros/pinv.sci
new file mode 100755
index 000000000..056ad2f66
--- /dev/null
+++ b/modules/linear_algebra/macros/pinv.sci
@@ -0,0 +1,42 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function X=pinv(A,tol)
+    //Pseudo inverse
+    if type(A)==1 then
+        if A==[] then X=[],return,end
+
+        [U,S,V] = svd(A,"e");
+        S = diag(S)
+        if argn(2) < 2
+            tol = max(size(A)) * S(1) * %eps;
+        end
+        r=size(find(S>tol),"*") //Rank
+        if r == 0 then
+            X = zeros(A');
+        else
+            //X=V1*inv(S1)*U1'
+            X = V(:,1:r)* diag(ones(r,1)./S(1:r)) *U(:,1:r)';
+        end
+    else
+        [t,n]=typename();n=stripblanks(n(find(t==type(A))))
+        fun="%"+n+"_pinv"
+        if exists(fun)==1 then
+            if argn(2)==1 then
+                execstr("X="+fun+"(A)")
+            else
+                execstr("X="+fun+"(A,tol)")
+            end
+        else
+            error(msprintf(gettext("%s: Function not defined for type ''%s''. Check argument or define function %s."),"pinv",n,fun));
+        end
+    end
+endfunction
+
diff --git a/modules/linear_algebra/macros/polar.bin b/modules/linear_algebra/macros/polar.bin
new file mode 100755
index 000000000..228b528e1
--- /dev/null
+++ b/modules/linear_algebra/macros/polar.bin
diff --git a/modules/linear_algebra/macros/polar.sci b/modules/linear_algebra/macros/polar.sci
new file mode 100755
index 000000000..0d416625a
--- /dev/null
+++ b/modules/linear_algebra/macros/polar.sci
@@ -0,0 +1,29 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [ro,teta]=polar(a);
+    //[ro,Teta]=polar(A);
+    // Polar form of A
+    // A=Ro*exp(%i*Teta) Ro symmetric >=0
+    // Teta hermitian >=0
+    //F.D.
+    //!
+    [U,s,V]=svd(a);
+    ro1=U*sqrt(s);
+    ro=ro1*ro1';
+    W=U*V';
+    // A = Ro*W   (Ro sdp ; W unit)
+    // W=exp(%i*Teta)
+    //
+    [ab,x,bs]=bdiag(W+0*%i*W);
+    z=log(diag(ab));
+    lw=x*diag(z)*inv(x);
+    teta=-%i*lw;
+endfunction
diff --git a/modules/linear_algebra/macros/proj.bin b/modules/linear_algebra/macros/proj.bin
new file mode 100755
index 000000000..aff7c53d2
--- /dev/null
+++ b/modules/linear_algebra/macros/proj.bin
diff --git a/modules/linear_algebra/macros/proj.sci b/modules/linear_algebra/macros/proj.sci
new file mode 100755
index 000000000..ccb93f635
--- /dev/null
+++ b/modules/linear_algebra/macros/proj.sci
@@ -0,0 +1,19 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [y]=proj(x1,x2)
+    //[y]=proj(x1,x2)  projection on x2 parallel to x1
+    //F.D.
+    //!
+    [l,k]=size(x1);
+    [w,n]=rowcomp(x1);w1=w(:,1:n);
+    x1t=w(n+1:l,:);
+    y=x2/(x1t*x2)*x1t
+endfunction
diff --git a/modules/linear_algebra/macros/projspec.bin b/modules/linear_algebra/macros/projspec.bin
new file mode 100755
index 000000000..424f4573e
--- /dev/null
+++ b/modules/linear_algebra/macros/projspec.bin
diff --git a/modules/linear_algebra/macros/projspec.sci b/modules/linear_algebra/macros/projspec.sci
new file mode 100755
index 000000000..f3080787c
--- /dev/null
+++ b/modules/linear_algebra/macros/projspec.sci
@@ -0,0 +1,61 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [S,P,D,index]=projspec(A,tol)
+    //[S,P,D,index]=projspec(A)
+    //Spectral characteristics of A at 0
+    //S = reduced resolvent at 0 (S=-Drazin_inverse(A))
+    //P = spectral projection at 0
+    //D = Nilpotent operator at 0
+    //index = index of the 0 eigenvalue
+    //!
+    [LHS,RHS]=argn(0)
+    [n,n]=size(A);
+    if RHS==1 then tol=1.d-6;end;
+    // A=0 ?
+    if norm(A,1) < %eps*n*n
+        P=eye(A),D=A,S=0*P;index=1;
+    end
+    //     nonsingular A:  index 0
+    if rcond(A) > tol
+        S=inv(A),P=0*eye(A);D=P;index=0;return;
+    end;
+    //write(%io(2),'    rank A^k    rcond')
+    //                    index 1
+    index=1;
+    [B,C,dim]=fullrf(A);
+    if dim==0
+        P=eye(n,n);S=0*P;D=A;return
+    end;
+    Ck=C;Bk=B;//write(%io(2),[dim,rcond(C*B)],'(7x,f3.0,6x,e9.3)');
+    tst=rcond(Ck*Bk);
+    if size(Ck,1)==1 then tst=norm(Ck*Bk,1);end
+    if tst > tol then
+        M=inv(C*B);P=eye(A)-B*M*C;S=B*M*M*C;D=0*A;return
+    end
+    //                  Higher index
+    for k=2:n
+        [B,C,dim]=fullrf(C*B);
+        if dim==0
+            P=eye(n,n);S=0*P;D=A;index=k;return;
+        end;
+        Bk=Bk*B;Ck=C*Ck;      // Bk*Ck = A^k  (Full rank factorization)
+        index=k;
+        //write(%io(2),[dim,rcond(C*B)],'(7x,f3.0,6x,e9.3)');
+        if norm(C*B)*rcond(C*B) > tol,
+            M=inv((C*B)**index);   //  M=inv(Ck*Bk);    (Alternative computation)
+            P=eye(n,n)-Bk*M*Ck;   // P=eye(n,n)-Bk*inv(Ck*Bk)*Ck;
+            S=Bk*M*inv(C*B)*Ck;   // S=inv(A-P-D)+P
+            D=0.5*(A*P+P*A);return;
+        end
+    end;
+    P=eye(n,n);S=0*P;D=A;index=k;return;
+endfunction
+
diff --git a/modules/linear_algebra/macros/psmall.bin b/modules/linear_algebra/macros/psmall.bin
new file mode 100755
index 000000000..a48865079
--- /dev/null
+++ b/modules/linear_algebra/macros/psmall.bin
diff --git a/modules/linear_algebra/macros/psmall.sci b/modules/linear_algebra/macros/psmall.sci
new file mode 100755
index 000000000..1002dee0f
--- /dev/null
+++ b/modules/linear_algebra/macros/psmall.sci
@@ -0,0 +1,50 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Q,M]=psmall(A,thres,flag)
+    // [Q,M]=psmall(A,thres,flag)
+    //Projection on eigensubspace associated with eigenvalues
+    //with real part < thres (case flag='c') or with modulus <thres (flag='d')
+    //Projection is defined by Q*M. Eigenvalues of M*A*Q = eigenvalues
+    //of A with real part < thres (case flag='c',...).
+    //If [Q1,M1]== full rank factorization (fullrf) of eye-Q*M then evals of
+    // M1*A*Q1 =evals of A with real part >= thres (case flag='d',...).
+    // See also pbig
+    //!
+    [n,n]=size(A);
+    thres=real(thres);
+    if flag=="c" then
+        deff("[flag]=%smallei(x)","flag=real(x) < thres")
+        deff("[flag]=%bigeig(x)","flag=real(x) >= thres")
+
+    elseif flag=="d" then
+        deff("[flag]=%smallei(x)","flag=abs(x) < thres")
+        deff("[flag]=%bigeig(x)","flag=abs(x) >= thres")
+    else
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"psmall",3,"''c'',''d''"));
+    end
+    //
+    [X,dsmall] = schur(A,%smallei);
+    [Y,dbig] = schur(A,%bigeig);
+    Q=X(:,1:dsmall);
+    if Q==[] then M=[];return;end
+    Y1=Y';
+    M1=Y1(dbig+1:n,:);
+    E=M1*Q;
+    if rcond(E)>1.d-6 then
+        M=E\M1;
+    else
+        [Ab,X0]=balanc(A);
+        [X,dsmall] = schur(Ab,%smallei);X1=X*X0;Q=X1(:,1:dsmall);
+        [Y,dbig] = schur(Ab,%bigeig);Y1=inv(X0)*Y';M=Y1(dbig+1:n,:);
+        E=M*Q;
+        M=E\M;
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/quaskro.bin b/modules/linear_algebra/macros/quaskro.bin
new file mode 100755
index 000000000..7419b9a37
--- /dev/null
+++ b/modules/linear_algebra/macros/quaskro.bin
diff --git a/modules/linear_algebra/macros/quaskro.sci b/modules/linear_algebra/macros/quaskro.sci
new file mode 100755
index 000000000..8a7417c7c
--- /dev/null
+++ b/modules/linear_algebra/macros/quaskro.sci
@@ -0,0 +1,52 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+
+function [Q,Z,Ec,Ac,Qd,Zd,numbeps]=quaskro(E,A,tol)
+    // quasi- Kronecker form: s*Ec - Ac = Q*(sE-A)*Z
+    //
+    //             | sE(eps)-A(eps) |        X       |      X     |
+    //             |----------------|----------------|------------|
+    //             |        O       | sE(inf)-A(inf) |      X     |
+    //  Q(sE-A)Z = |=================================|============|
+    //             |                                 |            |
+    //             |                O                | sE(r)-A(r) |
+    //
+    // Ec=Q*E*Z, Ac=Q*A*Z, eps=Qd(1) x Zd(1) ,inf=Qd(2) x Zd(2)
+    // r = Qd(3) x Zd(3)
+    // numbeps(1) = # of eps blocks of size 0 x 1
+    // numbeps(2) = # of eps blocks of size 1 x 2
+    // numbeps(3) = # of eps blocks of size 2 x 3     etc...
+    // interface  from Slicot-fstair (F.D.)
+    // T. Beelen's routines
+    //!
+    [LHS,RHS]=argn(0);
+    if RHS==1 then [E,A]=pen2ea(E);tol=1.d-10;end
+    if RHS==2 then
+        if type(E)==2 then [E,A]=pen2ea(E);end  //quaskro(pencil,tol)
+        if type(E)==1 then tol=1.d-10;end   //quaskro(E,A);
+    end
+    [na,ma]=size(A);
+    Q=eye(na,na);Z=eye(ma,ma);
+    if E~=[] then nE=norm(E,1);else nE=0;end
+    [E,Q,Z,stair,rk]=ereduc(E,1000*%eps+tol*nE)
+    A=Q*A*Z;
+
+    if A~=[] then
+        tol=tol*max([norm(A,"fro"),norm(E,"fro")])+10*tol;
+    else
+        tol=0
+    end
+    [Ac,Ec,Q,Z,nlbcks,muk,nuk,muk0,nuk0,mnei]=fstair(A,E,Q,Z,stair,rk,tol)
+    numbeps=muk0(1:nlbcks)-nuk0(1:nlbcks);
+    Qd=[mnei(1),mnei(3),na-mnei(1)-mnei(3)];
+    Zd=[mnei(2),mnei(3),ma-mnei(2)-mnei(3)];
+
+endfunction
diff --git a/modules/linear_algebra/macros/randpencil.bin b/modules/linear_algebra/macros/randpencil.bin
new file mode 100755
index 000000000..08d3d1a4b
--- /dev/null
+++ b/modules/linear_algebra/macros/randpencil.bin
diff --git a/modules/linear_algebra/macros/randpencil.sci b/modules/linear_algebra/macros/randpencil.sci
new file mode 100755
index 000000000..28d01d1c2
--- /dev/null
+++ b/modules/linear_algebra/macros/randpencil.sci
@@ -0,0 +1,95 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function F=randpencil(eps,infi,fin,eta)
+    //returns a random pencil with given kronecker structure
+    //eps=[eps1,...,epsk]; epsilon blocks (size eps1x(eps1+1),....)
+    //fin=[l1,...,ln]  finite eigenvalues (assumed real)  (possibly [])
+    //infi=[k1,...,kp] size of J blocks, ki>=1  (infi=[] if no J blocks)
+    //eta=[eta1,...,etap];   eta blocks size (eta1+1)xeta1,...)
+    // epsi's should be >=0
+    // etai's should be >=0
+    // infi's should be >=1
+    // If called with eps=[0,...,0], infi=[], fin=[], eta=[]
+    // randpencil returns F=[];
+    // this should be an empty matrix with zero rows and coldim(eps) columns
+    // If called with eps=[], infi=[], fin=[], eta=[0,..,0]
+    // randpencil returns F=[];
+    // this should be an empty matrix with coldim(eta) rows and 0 columns.
+    // (bad behavior of the empty matrix!!!!!)
+
+    [LHS,RHS]=argn(0);
+    if RHS<>4 then
+        error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),"randpencil",4));
+    end
+    select type(fin)
+    case 1
+        if diag(fin)<>[] then Fin=%s*eye()-diag(fin);else Fin=[];end
+    case 2
+        Fin=%s*eye()-companion(fin);
+    end
+    deff("j=%jdrn(n)","j=zeros(n,n);for k=1:n-1;j(k,k+1)=1;end")
+    deff("Lk=%eta(k)","Lk=zeros(k+1,k);if k==0 then return;end;for j=1:k;Lk(j,j)=%s;Lk(j+1,j)=-1;end");
+    deff("Lk=%epsilon(k)","Lk=zeros(k,k+1);if k==0 then return;end;for j=1:k;Lk(j,j)=%s;Lk(j,j+1)=-1;end");
+
+    J=[];
+    for kk=infi;
+        J=sysdiag(J,%jdrn(kk));
+    end
+    if J==[] then Infin=[],else Infin=%s*J-eye();end
+
+    flageps=%f;
+    Eps=[];
+    seps=gsort(eps);
+    if seps(1)==0 then flageps=%t;end
+    if ~flageps then
+        for k=seps;
+            if k==0 then [p,q]=size(Eps); Eps=[Eps,zeros(p,1)];end
+            if k<>0 then Eps=sysdiag(Eps,%epsilon(k));end
+        end
+    end
+
+    flageta=%f;
+    Eta=[];
+    seta=gsort(eta);
+    if seta(1)==0 then flageta=%t;end
+    if ~flageta then
+        for k=seta;
+            if k==0 then [p,q]=size(Eta); Eta=[Eta;zeros(1,q)];end
+            if k<>0 then Eta=sysdiag(Eta,%eta(k));end
+        end
+    end
+
+    F=sysdiag(Eps,Infin,Fin,Eta);
+
+    [p,q]=size(F);ncols=q;
+
+    if flageps then
+        F=[zeros(p,prod(size(eps))),F];
+        if F==[] then ncols=prod(size(eps));end
+    end
+
+    if flageta then
+        [p,q]=size(F);
+        if F~=[] then
+            F=[F;zeros(prod(size(eta)),q)];
+        else
+            F=[F;zeros(prod(size(eta)),ncols)];
+        end
+    end
+    // This can be uncommented for a seemingly more random pencil!
+    //[p,q]=size(F);
+    //rand('seed',0);
+    //rand('normal')
+    //Q=rand(p,p);
+    //Z=rand(q,q);
+    //F=Q*F*Z;
+endfunction
+
diff --git a/modules/linear_algebra/macros/range.bin b/modules/linear_algebra/macros/range.bin
new file mode 100755
index 000000000..69ce81568
--- /dev/null
+++ b/modules/linear_algebra/macros/range.bin
diff --git a/modules/linear_algebra/macros/range.sci b/modules/linear_algebra/macros/range.sci
new file mode 100755
index 000000000..ed45a765b
--- /dev/null
+++ b/modules/linear_algebra/macros/range.sci
@@ -0,0 +1,34 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [X,dim]=range(A,k)
+    // Computation of Range A^k ; the first dim rows of X span the
+    // range of A^k.
+    //!
+    if argn(2)==1 then k=1,end
+    k=double(k)
+    if int(k)<>k|k<0 then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Non-negative integer expected.\n"),"range",2));
+    end
+    if size(A,1)<>size(A,2)|~isreal(A) then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"range",1));
+    end
+
+    if k==0 then
+        dim=size(A,1);X=eye(A);
+    else
+        [U,dim]=rowcomp(A);X=U;
+        for l=2:k
+            A=A*U';
+            [U,dim]=rowcomp(A(:,1:dim));
+            X=U*X;
+        end;
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/rank.bin b/modules/linear_algebra/macros/rank.bin
new file mode 100755
index 000000000..a0c915c24
--- /dev/null
+++ b/modules/linear_algebra/macros/rank.bin
diff --git a/modules/linear_algebra/macros/rank.sci b/modules/linear_algebra/macros/rank.sci
new file mode 100755
index 000000000..9de89b22c
--- /dev/null
+++ b/modules/linear_algebra/macros/rank.sci
@@ -0,0 +1,33 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function r=rank(A,tol)
+    //Matrix rank
+    if type(A)==1 then
+        if A==[] then r=0;return,end
+        s = svd(A);
+        if argn(2) == 1 then
+            tol = max(size(A)) * s(1) * %eps;
+        end
+        r = size(find(s > tol),"*");
+    else
+        [t,n]=typename();n=stripblanks(n(find(t==type(A))))
+        fun="%"+n+"_rank"
+        if exists(fun)==1 then
+            if argn(2)==1 then
+                execstr("r="+fun+"(A)")
+            else
+                execstr("r="+fun+"(A,tol)")
+            end
+        else
+            error(msprintf(gettext("%s: Function not defined for type ''%s''. Check argument or define function %s."),"rank",n,fun));
+        end
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/rowcomp.bin b/modules/linear_algebra/macros/rowcomp.bin
new file mode 100755
index 000000000..42e1d1578
--- /dev/null
+++ b/modules/linear_algebra/macros/rowcomp.bin
diff --git a/modules/linear_algebra/macros/rowcomp.sci b/modules/linear_algebra/macros/rowcomp.sci
new file mode 100755
index 000000000..d3292f263
--- /dev/null
+++ b/modules/linear_algebra/macros/rowcomp.sci
@@ -0,0 +1,46 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 1987-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [w,rk]=rowcomp(A,flag,tol)
+    //Row compression of A <--> computation of im(A)
+    //flag and tol are optional parameters
+    //flag='qr' or 'svd' (default 'svd')
+    //tol tolerance parameter (sqrt(%eps)*norm(A,1) as default value)
+    //the rk first (top) rows of w span the row range of a
+    //the rk first columns of w' span the image of a
+
+    if A==[] then w=[];rk=0;return;end
+
+    [ma,na]=size(A)
+    rhs=argn(2)
+
+    if norm(A,1) < sqrt(%eps)/10 then rk=0,w=eye(ma,ma),return;end
+
+    select rhs
+    case 1 then// default values for flag and tol
+        flag="svd",tol=sqrt(%eps)*norm(A,1);
+    case 2 then //default value for tol
+        tol=sqrt(%eps)*norm(A,1)
+    else
+        if size(tol,"*")>1|~isreal(tol)|tol<0 then
+            error(msprintf(gettext("%s: Wrong values for input argument #%d: Non-negative scalar expected.\n"),"rowcomp",3));
+        end
+    end
+    select flag
+    case "qr" then
+        [q,r,rk,e]=qr(A,tol);w=q';
+    case "svd" then
+        [u,s,v,rk]=svd(A,tol);w=u' ;
+    else
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"rowcomp",2,"''qr'',''svd''"));
+    end
+endfunction
+
+
diff --git a/modules/linear_algebra/macros/rowshuff.bin b/modules/linear_algebra/macros/rowshuff.bin
new file mode 100755
index 000000000..39adf82a6
--- /dev/null
+++ b/modules/linear_algebra/macros/rowshuff.bin
diff --git a/modules/linear_algebra/macros/rowshuff.sci b/modules/linear_algebra/macros/rowshuff.sci
new file mode 100755
index 000000000..42dabd9fd
--- /dev/null
+++ b/modules/linear_algebra/macros/rowshuff.sci
@@ -0,0 +1,50 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Ws,Fs1]=rowshuff(Fs,alfa)
+    // Shuffle algorithm: Given the pencil Fs=s*E-A, returns Ws=W(s)
+    // (square polynomial matrix) such that:
+    // Fs1 = s*E1-A1 = W(s)*(s*E-A) is a pencil with non singular E1 matrix.
+    // This is possible iff the pencil Fs = s*E-A is regular (i.e. invertible).
+    // The poles @ infinity of Fs are put to alfa and the zeros of Ws are @ alfa.
+    // Note that (s*E-A)^-1 = (s*E1-A1)^-1 * W(s) = (W(s)*(s*E-A))^-1 *W(s)
+
+    [LHS,RHS]=argn(0);
+    if RHS==1 then
+        alfa=0;
+    end
+    [E,A]=pen2ea(Fs);
+    //     E is non singular: --> exit
+    if rcond(E) >= 1.d-5 then W=eye(E);Fs1=Fs;return;end
+    //     E is singular:
+    s=poly(0,"s");tol=1.d-10*(norm(E,1)+norm(A,1));
+    [n,n]=size(E);Ws=eye(n,n);
+    //
+    rk=0;i=0;
+    while rk  < n
+        if i==n then
+            error(msprintf(gettext("%s: Singular pencil."),"rowshuffle"));
+            W=[];
+        end
+        [W,rk]=rowcomp(E);
+        if rk==n then return;end
+        W1=W(1:rk,:);W2=W(rk+1:n,:);
+        E=[W1*E;
+        -W2*A];
+        A=[W1*A;
+        -alfa*W2*A];
+        Fs1=s*E-A;
+        //     Update
+        Ws=[eye(rk,rk),zeros(rk,n-rk);
+        zeros(n-rk,rk),(s-alfa)*eye(n-rk,n-rk)]*W*Ws;
+        i=i+1;
+    end
+endfunction
+
diff --git a/modules/linear_algebra/macros/rref.bin b/modules/linear_algebra/macros/rref.bin
new file mode 100755
index 000000000..efe0eea80
--- /dev/null
+++ b/modules/linear_algebra/macros/rref.bin
diff --git a/modules/linear_algebra/macros/rref.sci b/modules/linear_algebra/macros/rref.sci
new file mode 100755
index 000000000..dc9bef402
--- /dev/null
+++ b/modules/linear_algebra/macros/rref.sci
@@ -0,0 +1,45 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [A,jb]=rref(A,tol)
+    //R = rref(A) produces the reduced row echelon form of A.
+    if argn(2)<2 then
+        tol=2*%eps*norm(A,"inf")*max(size(A))
+    end
+
+    if type(A)==1 then
+        if A==[] then R=[],return,end
+        [m,n]=size(A)
+        k = 1;l = 1;jb = [];
+        while (k <= m) & (l <= n)
+            // Find value and index of largest element in the remainder of column l.
+            [p,i] = max(abs(A(k:$,l))); i = i+k-1;
+            if (p <= tol) // The column is negligible
+                A(k:$,l) = zeros(m-k+1,1);
+                l = l + 1;
+            else
+                jb = [jb, l]; // Remember column index
+                A([k,i],l:$) = A([i,k],l:$); //swap rows i and j
+                A(k,l:$) = A(k,l:$)/A(k,l); // Normalize the pivot row
+
+                i = [1:k-1 k+1:m]; //other rows
+                A(i,l:$) = A(i,l:$) - A(i,l)*A(k,l:$); //eliminate
+                k = k + 1;l = l + 1;
+            end
+        end
+    else
+        [t,n]=typename();n=stripblanks(n(find(t==type(A))))
+        fun="%"+n+"_rref"
+        if exists(fun)==1 then
+            execstr("[A,jb]="+fun+"(A)")
+        else
+            error(msprintf(gettext("%s: Function not defined for type ''%s''. Check argument or define function %s."),"rref",n,fun));
+        end
+    end
+endfunction
diff --git a/modules/linear_algebra/macros/spaninter.bin b/modules/linear_algebra/macros/spaninter.bin
new file mode 100755
index 000000000..81d2f4737
--- /dev/null
+++ b/modules/linear_algebra/macros/spaninter.bin
diff --git a/modules/linear_algebra/macros/spaninter.sci b/modules/linear_algebra/macros/spaninter.sci
new file mode 100755
index 000000000..e181e4a05
--- /dev/null
+++ b/modules/linear_algebra/macros/spaninter.sci
@@ -0,0 +1,52 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [x,dim]=spaninter(a,b,tol)
+    //[X,dim]=spaninter(a,b [,tol])  computes intersection of Range(A)
+    // and Range(B)
+    //
+    // x(:,1:dim) is an orthogonal basis for A inter B.
+    //            In the X basis A and B are respectively:
+    //            X'*A and X'*B.
+    // dim        dimension of subspace A inter B.
+    // tol        threshold (sqrt(%eps) is the default value).
+
+    [lhs,rhs]=argn(0);
+    [na,ma]=size(a);[nb,mb]=size(b);
+    if ma*na==0 then dim=0;x=eye(nb,nb);return;end
+    if mb*nb==0 then dim=0;x=eye(na,na);return;end
+    if rhs==2 then tol=sqrt(%eps);end
+    if na <> nb then
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"),"spaninter",1,2));
+    end
+    if mb > ma then [x,dim]=spaninter(b,a,tol),return,end
+    [xp,ra]=rowcomp(a);x=xp'
+    //test  trivial cases :
+    //a surjective a inter b is b
+    if ra == na then [xp,dim]=rowcomp(b),x=xp',return,end
+    //a is zero
+    if ra == 0 then dim=0,return,end
+    b=xp*b;      //update
+
+    // b2=vectors in b whch are in a
+    up=1:ra;low=ra+1:na;
+    [v1,k2]=colcomp(b(low,:));
+    b1=b*v1;     //update
+    bup=b1(up,:);blow=b1(low,:)
+    if norm(blow,1) <= tol*norm(bup,1) then k2=0,end
+    k1=mb-k2;
+    if k1==0 then dim=0;return;end
+    b2=b1(:,1:k1);  //intersection in base x
+    if norm(b2,1) < tol*norm(b,1)*mb*nb then dim=0,return,end
+    [u2p,dim]=rowcomp(b2(up,:));
+    x(:,up)=x(:,up)*u2p'         //update
+endfunction
+
+
+
diff --git a/modules/linear_algebra/macros/spanplus.bin b/modules/linear_algebra/macros/spanplus.bin
new file mode 100755
index 000000000..45e19018d
--- /dev/null
+++ b/modules/linear_algebra/macros/spanplus.bin
diff --git a/modules/linear_algebra/macros/spanplus.sci b/modules/linear_algebra/macros/spanplus.sci
new file mode 100755
index 000000000..2f2f33be6
--- /dev/null
+++ b/modules/linear_algebra/macros/spanplus.sci
@@ -0,0 +1,40 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [x,dim,dima]=spanplus(a,b,tol)
+    //[X,dim,dima]=spanplus(A,B,tol) computes an orthogonal basis of
+    // a+b such that : the first dima columns of x span Range(A)
+    // and the following (dim-dima) columns make a basis of a+b
+    // relative to a. tol is an optional argument.
+    // The dim first columns of x make a basis for A+B.
+    //!
+    [na,ma]=size(a);[nb,mb]=size(b);
+    if na*ma==0 then
+        dima=0;[x,dim]=rowcomp(b);x=x';return;
+    end
+    if nb*mb==0 then [x,dima]=rowcomp(a);dim=dima;x=x';return;end
+    [lhs,rhs]=argn(0);
+    if rhs==2 then tol=%eps*na*nb*ma*mb;end
+    if na<>nb then
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"),"spanplus",1,2));
+    end
+    [x,dima]=rowcomp(a);
+    b=x*b;x=x'    //update b,x
+    if dima == na then dim=dima,return,end;
+    low=(dima+1):na;
+    blow=b(low,:);
+    if norm(blow,1) <= tol*norm(b,1),dim=dima,return,end
+    [u2,r2]=rowcomp(blow);
+    dim=dima+r2;
+    x(:,low)=x(:,low)*u2';    //update
+endfunction
+
+
+
diff --git a/modules/linear_algebra/macros/spantwo.bin b/modules/linear_algebra/macros/spantwo.bin
new file mode 100755
index 000000000..b2b50afba
--- /dev/null
+++ b/modules/linear_algebra/macros/spantwo.bin
diff --git a/modules/linear_algebra/macros/spantwo.sci b/modules/linear_algebra/macros/spantwo.sci
new file mode 100755
index 000000000..f52a36efb
--- /dev/null
+++ b/modules/linear_algebra/macros/spantwo.sci
@@ -0,0 +1,48 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA - François DELEBECQUE
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [Xp,dima,dimb,dim]=spantwo(A,B)
+    //Given two matrices A and B with same number of rows,
+    //returns a square matrix Xp (not necessarily orthogonal)
+    //such that :
+    //         [*,0]    (dim-dimb rows)
+    //Xp*[A,B]=[*,*]    (dima+dimb-dim rows)
+    //         [0,*]    (dim-dima rows)
+    //         [0,0]
+    //The  dima first columns of inv(Xp) span range(A).
+    //Columns dim-dimb+1 to dima of inv(Xp) span the intesection of
+    //range(A) and range(B).
+    //The dim first columns of inv(Xp) span range(A)+range(B).
+    // Ex: A=[1,0,0,4;
+    //        5,6,7,8;
+    //        0,0,11,12;
+    //        0,0,0,16];
+    //B=[1,2,0,0]';C=[4,0,0,1];
+    //Sl=ss2ss(syslin('c',A,B,C),rand(A));
+    //[no,X]=contr(Sl(2),Sl(3));CO=X(:,1:no);  //Controllable part
+    //[uo,Y]=unobs(Sl(2),Sl(4));UO=Y(:,1:uo);  //Unobservable part
+    //[Xp,dimc,dimu,dim]=spantwo(CO,UO);    //Kalman decomposition
+    //Slcan=ss2ss(Sl,inv(Xp));
+
+    [X1,dim,dima]=spanplus(A,B);Xp=X1';
+    B1B2B3=Xp*B;B1B2B3=B1B2B3(1:dim,:);
+    [W,dimb]=rowcomp(B1B2B3);W=W(dim:-1:1,:);
+    W11=W(1:dima,1:dima);W21=W(dima+1:dim,1:dima);
+    if rcond(W11)<1.d-10 then
+        //   Which is better?
+        B1B2=B1B2B3(1:dima,:);[W,dimb0]=rowcomp(B1B2);W=W(dima:-1:1,:);
+        [n1,n2]=size(A);
+        Xp=sysdiag(W,eye(n1-dima,n1-dima))*Xp;
+        return;
+    end
+    Q=[eye(dima,dima),zeros(dima,dim-dima);
+    -W21*inv(W11),eye(dim-dima,dim-dima)];
+    Xp(1:dim,:)=Q*W*Xp(1:dim,:);
+endfunction
diff --git a/modules/linear_algebra/macros/sqroot.bin b/modules/linear_algebra/macros/sqroot.bin
new file mode 100755
index 000000000..c17947bd7
--- /dev/null
+++ b/modules/linear_algebra/macros/sqroot.bin
diff --git a/modules/linear_algebra/macros/sqroot.sci b/modules/linear_algebra/macros/sqroot.sci
new file mode 100755
index 000000000..88a9253cf
--- /dev/null
+++ b/modules/linear_algebra/macros/sqroot.sci
@@ -0,0 +1,24 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [S]=sqroot(Q)
+    Q1=(Q+Q')/2;
+    if norm(Q1-Q,1) > 100*%eps then
+        warning(msprintf(gettext("%s: Wrong size for input argument #%d: Symmetric expected"),"sqroot",1));
+    end
+    tt=min(spec(Q1));
+    if tt <-10*%eps then
+        warning(msprintf(gettext("%s: Wrong value for input argument #%d: Not semi-definite positive"),"sqroot",1));
+    end
+    if norm(Q,1) < sqrt(%eps) then S=[];return;end
+    [u,S,v,rk]=svd(Q);
+    S=v(:,1:rk)*sqrt(S(1:rk,1:rk));
+    if norm(imag(Q1),1) <1.d-8 then S=real(S);;end
+endfunction
diff --git a/modules/linear_algebra/macros/squeeze.bin b/modules/linear_algebra/macros/squeeze.bin
new file mode 100755
index 000000000..e4164b68a
--- /dev/null
+++ b/modules/linear_algebra/macros/squeeze.bin
diff --git a/modules/linear_algebra/macros/squeeze.sci b/modules/linear_algebra/macros/squeeze.sci
new file mode 100755
index 000000000..0a2114d9f
--- /dev/null
+++ b/modules/linear_algebra/macros/squeeze.sci
@@ -0,0 +1,37 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2002-2007 - INRIA - Eric Dubois
+// Copyright (C) 2007 - INRIA - Jean-Baptiste Silvy
+// Copyright (C) 2011 - INRIA - Serge Steer (extension to Structs, Cells,
+//                                           2D arrays of any types)
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function A = squeeze(A)
+
+    // PURPOSE: Remove singleton dimensions, that is any dimension
+    // for which the size of the input hypermatrix is 1; if the
+    // input is a matrix, it is unaffected
+    // ------------------------------------------------------------
+    // INPUT:
+    // * A = a hypermatrix or a matrix
+    // ------------------------------------------------------------
+
+    if or(typeof(A)==["hypermat","ce","st"]) then
+        Dims=size(A)
+        newDims = Dims(Dims <> 1) ;
+        if size(newDims,"*") <2  then
+            A=A(:)
+        else
+            A=matrix(A,newDims)
+        end
+    elseif type(A)<=10 then
+        // it is a standard matrix nothing to do
+    else
+        error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"squeeze",1))
+    end
+
+endfunction
diff --git a/modules/linear_algebra/macros/sva.bin b/modules/linear_algebra/macros/sva.bin
new file mode 100755
index 000000000..de527bf20
--- /dev/null
+++ b/modules/linear_algebra/macros/sva.bin
diff --git a/modules/linear_algebra/macros/sva.sci b/modules/linear_algebra/macros/sva.sci
new file mode 100755
index 000000000..98135b85a
--- /dev/null
+++ b/modules/linear_algebra/macros/sva.sci
@@ -0,0 +1,32 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [U,S,V]=sva(A,tol)
+    // approximation de rang donne d'une matrice
+    if A==[] then U=[],S=[],V=[],return,end
+    [U,S,V]=svd(A,"e")
+    if argn(2)==1 then
+        tol = max(size(A)) * S(1) * %eps;
+        rk = size(find(diag(S) > tol),"*");
+    else
+        if tol>1 then //rank given
+            rk=tol
+            if rk>min(size(A)) then
+                error(msprintf(gettext("%s: Wrong value for input argument #%d: Requested rank is greater than matrix dimension."),"sva",1));
+            end
+        else
+            rk = size(find(diag(S) > tol),"*");
+        end
+    end
+    U=U(:,1:rk);S=S(1:rk,1:rk),V=V(:,1:rk)
+endfunction
+
+
+
diff --git a/modules/linear_algebra/macros/trace.bin b/modules/linear_algebra/macros/trace.bin
new file mode 100755
index 000000000..1ae164f0f
--- /dev/null
+++ b/modules/linear_algebra/macros/trace.bin
diff --git a/modules/linear_algebra/macros/trace.sci b/modules/linear_algebra/macros/trace.sci
new file mode 100755
index 000000000..953502bc2
--- /dev/null
+++ b/modules/linear_algebra/macros/trace.sci
@@ -0,0 +1,56 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function t=trace(a)
+    // trace - computes the trace of a matrix
+    select type(a)
+    case 1 then
+        [m,n]=size(a)
+        if m<>n then
+            error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"trace",1));
+        end
+        t=sum(diag(a))
+    case 2 then
+        [m,n]=size(a)
+        if m<>n then
+            error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"trace",1));
+        end
+        t=sum(diag(a))
+        //-compat next case retained for list /tlist compatibility
+    case 15 then
+        if a(1)=="r" then
+            [m,n]=size(a)
+            if m<>n then
+                error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"trace",1));
+            end
+            t=sum(diag(a))
+        else
+            error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"trace",1));
+        end
+    case 16 then
+        if a(1)=="r" then
+            [m,n]=size(a)
+            if m<>n then
+                error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"trace",1));
+            end
+            t=sum(diag(a))
+        else
+            error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"trace",1));
+        end
+    case 5 then
+        [m,n]=size(a)
+        if m<>n then
+            error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"trace",1));
+        end
+        t=sum(diag(a))
+    else
+        error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"trace",1));
+    end
+endfunction