function [ok] = newton_step_L1_ok()
    ok = true;
    % --------------------------------------------------------
    % You can change these parameters
    m    = 2;   % number of measurements per time point
    n    = 3;   % number of state vector components per time point
    N    = 4;   % number of time points
    % ---------------------------------------------------------
    %  Define the problem
    rand('seed', 123);
    mu    = .5;
    s     = rand(m, N);
    y     = rand(n, N);
    r     = rand(m, N);
    b     = rand(m, N);
    d     = rand(n, N);
    Hdia  = zeros(n, n, N);
    Hlow  = zeros(n, n, N);
    Bdia  = zeros(m, n, N);
    pPlus = rand(m, N);
    pMinus = rand(m, N);
    for k = N : -1 : 1
        Bdia(:,:, k)    = rand(m, n);
        tmp             = rand(n, n);
        Hlow(:,:, k)    = tmp;
        Hdia(:,:, k)   = (tmp * tmp') + 4 * eye(n);
        %
    end

    % Form full matrices.
    H = ckbs_blktridiag_mul(Hdia, Hlow, speye(n*N));
    B = ckbs_blkdiag_mul(Bdia, speye(n*N));

    % -------------------------------------------------------------------
    [dPPlus, dPMinus, dR, dS, dY] = ckbs_newton_step_L1(mu, s, y, r, b, d, Bdia, Hdia, Hlow, pPlus, pMinus);
    % -------------------------------------------------------------------
    F      =ckbs_kuhn_tucker_L1(mu, s, y, r, b, d, Bdia, Hdia, Hlow, ...
        pPlus, pMinus);

    dF     = [ ...
        speye(m*N) , -speye(m*N), zeros(m*N), zeros(m*N), -B
        zeros(m*N) , diag(s(:)) , zeros(m*N), diag(pMinus(:)),zeros(m*N,n*N)
        zeros(m*N) , zeros(m*N) , speye(m*N), speye(m*N),zeros(m*N,n*N)
        diag(r(:)) , zeros(m*N) , diag(pPlus(:)), zeros(m*N),zeros(m*N,n*N)
        zeros(n*N, m*N), zeros(n*N, m*N), -B'/2, B'/2, H
        ];
    delta        = [ dPPlus(:); dPMinus(:); dR(:); dS(:); dY(:) ];
    ok           = ok & ( max( abs( F + dF * delta ) ) <= 1e-10 );
    return
end