function V = generate_V(Vchoice,n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% V is a basis matrix for the null space of ones(n,1)
% this function generated V in three different ways

if Vchoice==1,
    fprintf('using Vchoice = 1 standard V with qr for Vhat \n')
    % using qr
    Veye = [eye(n-1); -ones(1,n-1)];
    [V,~] = qr(Veye,0);
elseif Vchoice==2,
    fprintf('using  Vchoice = 2 uppertriangular V for Vhat \n')
    % using upper triangular ones
    V0 = ones(n,n-1);
    V0 = triu(V0);
    V0(2:n+1:end) = -(1:n-1);
    DV0 = sqrt((1:1:(n-1))+(1:1:(n-1)).^2);
    V0 = V0/diag(DV0);
    V = V0;
elseif Vchoice==3,
    fprintf('using  Vchoice = 3 sparsest V by blocks for Vhat \n')
    % using orthog. blocks
    iblk=1;                            % number of current block
    ii=1;                              % current beginning column of block iblk
    temp=(spalloc(n,n-1,n*(n+1)/2));     % matrix for final V
    sizeblocks=zeros(floor(n/2),1);    % initialize vector
    normalize=zeros(floor(n/2),1);    % initialize vector
    dscale=ones(n-1,1);              % scaling vector
    sizeblocks(iblk)=floor(n/(2^iblk));     % size/#cols of blk iblk
    normalize(iblk)=norm(ones(2^iblk,1));   %  normalization scalars
    while sizeblocks(iblk) >= 1,         % room for at least one column
        tspeye=speye(sizeblocks(iblk));
        ncols=sizeblocks(iblk);    
        % do not exceed # cols left % TK: this will never exceed cols left
        dscale(ii:ii+ncols-1)=normalize(iblk);

        temp(:,ii:ii+ncols-1)=...
            [ kron(tspeye(:,1:ncols),[ones(2^(iblk-1),1);-ones(2^(iblk-1),1)])
            zeros(n-sizeblocks(iblk)*2^iblk,ncols) ];
        ii=ii+ncols;
        iblk=iblk+1;
        sizeblocks(iblk)=floor(n/(2^iblk));
        normalize(iblk)=norm(ones(2^iblk,1));
    end
    itest= (sum(diag(temp'*temp)==0));
    if itest > 0,
        ttemp=sparse(null(full([temp ones(n,1)]')));
        temp(:,end-itest+1:end)=ttemp(:,1:itest);
    end
    dd=sqrt(diag(temp'*temp));
    Ddd=spdiags(dd,0,n-1,n-1);
    V=temp/Ddd;
end    % end if for Vchoice
if issparse(V)
    V = full(V);
end