function [ Q, fnew, iter, time] = NPG_constant(B, Q_0, L_0, s, opts)
% This function solves the following optimization problem by proximal gradient with majorization,
% min 1/2d^2(BQ,D),
% s.t. Q belong to C, C is a set of all r by r orthogonal matrix
% D = {B : max_{j=1,..,n} ||B_j||_0 <= s}, where B_j is the j-th row of B, B is an n by n matrix

% Input 
% B - n by n matrix 
% Q_0 - r by r orthogonal matrix
% L_0 - a number greater than lambda_max(B'B)
% s - sparsity level
% opts - contains the following structure variables, default values are in brackets:
%         
%           freq    - frequency for displaying output [5000]
%           maxiter - maximum number of iterations [inf]
%           disp    - display [1]: 1 - display per freq iterations
%                                  0 - no display
%           M       - memory in nonmonotone linesearch [5]
%

% Output
% 
% Q - approximate solution at termination
% iter - number of iterations
% fnew - function value at termination
% time - CPU time

tstart = tic;

% Parameters
if isfield(opts,'maxiter'), maxiter = opts.maxiter; else maxiter = 5000; end %5000, otherwise declare failure
if isfield(opts,'freq'), freq = opts.freq; else freq =5000; end
if isfield(opts,'disp'), disp = opts.disp; else disp = 1; end
if isfield(opts,'tol'), tol = opts.tol; else tol = 1e-20; end

% Set matrix handles
Ahandle = @(x) B*x;
AThandle = @(x) B'*x;

% Initialization
Q = Q_0;
iter = 1;

% Compute function value
tmp = Ahandle(Q);
tmp_1 = tmp - sparsematrix(tmp,s);

% compute gradient
g = AThandle(tmp_1);

%% Main interation
while 1 == 1
    
    z_2 = Q - 1/L_0*g;
    [U,~,V] = svd(z_2,'econ');
    Q = U*V';
    tmp = Ahandle(Q);
    tmp_1 = tmp - sparsematrix(tmp,s) ;
    fnew = 0.5*norm(tmp_1,'fro')^2;
    
    %% Display
    if mod(iter,freq) == 0 && disp == 1
        fprintf('iter = %5d, fval = %5.4e\n', iter, fnew);
    end  
    
    %% Check for termination
    if iter > maxiter
        fprintf('NPG fails\n')
        break
    end
    if sqrt(fnew*2) < tol
        fprintf('Successful termination\n')
        break
    end
    
    %% Update gradient and iterate
    g = AThandle(tmp_1);  
    iter = iter + 1;
end
time = toc(tstart);
fprintf('iter = %d, fnew = %g, time = %g\n', iter, fnew, time)
    
    