function [x0,V,rcnd]=inistate(A,B,C,D,y,u,tol,printw)
    x0=[];V=[];rcnd=[];
    [nargout,nargin] = argn(0)
    //INISTATE  Estimates the initial state of a discrete-time system, given the
    //          (estimated) system matrices, and a set of input/output data.
    //
    //        X0 = INISTATE(SYS,Y,U,TOL,PRINTW)  estimates the initial state X0 of
    //        the discrete-time system SYS = (A,B,C,D), using the output data Y
    //        and the input data U. The model structure is :
    //
    //             x(k+1) = Ax(k) + Bu(k),   k >= 1,
    //             y(k)   = Cx(k) + Du(k),
    //
    //        The vectors y(k) and u(k) are transposes of the k-th rows of Y and U,
    //        respectively.
    //        Instead of the first input parameter SYS (an ss object), equivalent
    //        information may be specified using matrix parameters, for instance,
    //        X0 = INISTATE(A,B,C,Y,U);   or   X0 = INISTATE(A,C,Y);
    //
    //        TOL is the tolerance used for estimating the rank of matrices.
    //        If  TOL > 0,  then the given value of  TOL  is used as a lower bound
    //        for the reciprocal condition number.
    //        Default:    prod(size(matrix))*epsilon_machine where epsilon_machine
    //                    is the relative machine precision.
    //
    //        PRINTW is a select for printing the warning messages.
    //        PRINTW = 1: print warning messages;
    //               = 0: do not print warning messages.
    //        Default:    PRINTW = 0.
    //
    //        [x0,V,rcnd] = INISTATE(SYS,Y,U,TOL,PRINTW)  returns, besides x0,
    //        the orthogonal matrix V which reduces the system state matrix A to
    //        a real Schur form, as well as an estimate of the reciprocal condition
    //        number of the coefficient matrix of the least squares problem solved.
    //
    //        See also FINDBD, FINDX0BD
    //

    //        V. Sima 13-05-2000.
    //
    //        For efficiency, most errors are checked in the mexfile findBD. Also,
    //        except for scalars, the input parameters are not copied, but renamed.
    //
    //        Revisions:
    //        V. Sima, July 2000.
    //

    ni = nargin;
    if mtlb_isa(A,"lti") then
        // Get the system matrices of the ss object, and the remaining parameters.
        // General call     x0 = inistate(A,B,C,D,y,u,tol,printw);
        // Special call     x0 = inistate(sys,y,u,tol,printw);
        //
        if A.dt=="c" then
            error(msprintf(gettext("%s: Wrong values for input argument #%d: Discrete time system expected.\n"),"inistate",1))
        end
        if ni<2 then
            error(msprintf(gettext("%s: Wrong number of input arguments: At least %d expected.\n"),"inistate",2));
        end
        [As,Bs,Cs,Ds] = abcd(A)
        [ny,nu] = size(A);
        [ty,p] = size(B);
        if ni>2 then
            [tu,m] = size(C);
            if ~(((tu==ty|tu==0)&(m==nu))&(ty>1)) then
                tol = C;
                // Special call     x0 = inistate(sys,y,tol,printw);
                if ni>3 then
                    printw = D;
                else
                    printw = 0;
                end
                // Below, B means y !
                [x0,Vl,rcndl] = findBD(1,3,As,Cs,B,tol,printw);
            else
                if ni>3 then
                    // Special call     x0 = inistate(sys,y,u,tol,printw);
                    tol = D;
                    if ni>4 then
                        printw = y;
                    else
                        printw = 0;
                    end
                else
                    tol = 0;
                    printw = 0;
                end
                // Below, B means y, and C means u !
                if norm(Ds,1)==0 then
                    [x0,Vl,rcndl] = findBD(1,2,1,As,Bs,Cs,B,C,tol,printw);
                else
                    [x0,Vl,rcndl] = findBD(1,2,2,As,Bs,Cs,Ds,B,C,tol,printw);
                end
            end
        else
            // Special call     x0 = inistate(sys,y);
            // Below, B means y !
            [x0,Vl,rcndl] = findBD(1,3,As,Cs,B);
        end
        //
    else
        // The system matrices are directly specified.
        // General call     x0 = inistate(A,B,C,D,y,u,tol,printw);
        // Special calls    x0 = inistate(A,B,C,y,u,tol,printw);
        //                  x0 = inistate(A,C,y,tol,printw);
        //
        if ni<3 then
            error(msprintf(gettext("%s: Wrong number of input arguments: At least %d expected.\n"),"inistate",3));
        end
        [m2,n2] = size(B);
        [m3,n3] = size(C);
        if ni>=4 then
            [m4,n4] = size(D);
            if ni>=5 then
                [m5,n5] = size(y);
                if ni>=6 then
                    [m6,n6] = size(u);
                    if ni>=7 then
                        if (ni==7)&(m6*n6>1) then
                            // Special call     x0 = inistate(A,B,C,D,y,u,tol);
                            [x0,Vl,rcndl] = findBD(1,2,1,A,B,C,D,y,u,tol);
                        elseif ni==7 then
                            // Special call     x0 = inistate(A,B,C,y,u,tol,printw);
                            // Below, D means y and y means u !
                            printw = tol;
                            tol = u;
                            [x0,Vl,rcndl] = findBD(1,2,1,A,B,C,D,y,tol,printw);
                        else
                            [x0,Vl,rcndl] = findBD(1,2,2,A,B,C,D,y,u,tol,printw);
                        end
                    else
                        if m6*n6>1 then
                            [x0,Vl,rcndl] = findBD(1,2,2,A,B,C,D,y,u);
                        else
                            // Special call     x0 = inistate(A,B,C,y,u,tol);
                            // Below, y means U and D means y !
                            tol = u;
                            [x0,Vl,rcndl] = findBD(1,2,1,A,B,C,D,y,tol);
                        end
                    end
                else
                    // Special calls    x0 = inistate(A,B,C,y,u);
                    //                  x0 = inistate(A,C,y,tol,printw);
                    if m5*n5>1 then
                        // Below, y means u and D means y !
                        [x0,Vl,rcndl] = findBD(1,2,1,A,B,C,D,y);
                    else
                        // Below, C means y and B means C !
                        tol = D;
                        printw = y;
                        [x0,Vl,rcndl] = findBD(1,3,A,B,C,tol,printw);
                    end
                end
            else
                // Below, D means tol, C means y, and B means C !
                [x0,Vl,rcndl] = findBD(1,3,A,B,C,D);
            end
        else
            // Below, C means y, and B means C !
            [x0,Vl,rcndl] = findBD(1,3,A,B,C);
        end
    end
    //
    if nargout>2 then
        V = Vl;
        rcnd = rcndl;
    elseif nargout>1 then
        V = Vl;
    end
    //
    // end inistate
endfunction