function [rho,Out] = GNQKD(MGamma,gamma,VKVlist,VZKVlist,opt)
%function [rho,Out] = GNQKD(MGamma,gamma,VKVlist,VZKVlist,opt)
%%% this code solves the (QKD) using GN method (after perhaps FR and rotation)
%%% prints out the optimality condition
%%%  (QKD)     min    f(rho) = VGV rho log VGV rho - VZKV rho log VZKV rho  (quantum relative entropy)
%%%            s.t.   MGamma rho(:) = gamma
%%%                   rho >= 0
%%% Input:   MGamma : vectorized constraint for Gamma(rho) = gamma
%%%          gamma : RHS vector for Gamma(rho) = gamma
%%%          VKVlist : cell array for kraus operator in VGV(rho)  
%%%          VZKVlist : cell array for  VZKV rho log VZKV rho
%%% Output: rho : optimal primal variable
%%%%        Out : various outputs containing dual variables, optimal value,
%%%%        <Notable fields in Out>
%%%%         Out.y = y: dual variable related to MGamma rho(:) = gamma
%%%%         Out.Z = Z: dual variable related to rho
%%%%         Out.bestub: best upper bound
%%%%         Out.bestlb: best lower bound
%%%%         Out.noiter: number of iterations

%% Copyrights
% Author: Hao Hu, Haesol Im and Henry Wolkowicz 
% Email: h92hu@uwaterloo.ca, j5im@uwaterloo.ca, hwolkowi@uwaterloo.ca
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% latest versionDec23.2020
%%%% latest versionJan12.2021
%%%% latest versionJan29.2021
%%%% latest versionMar19.2021
%%%% latest versionMar30.2021
%%%% latest versionJun5.2021
%%%% latest versionJun23.2021
%%%% latest versionJul30.2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if isfield(opt,'profilef'), profilef = opt.profilef;  else, profilef = 'noprofile';  end
if strcmp(profilef,'solver')
    profile clear
    profile on -timestamp
end

if isfield(opt,'tolerGN'), tolerGN = opt.tolerGN;  else, tolerGN = 1e-12;  end
%% crosstoler for checking upper and lower bounds availability
if isfield(opt,'crosstoler'), crosstoler = opt.crosstoler;  else, crosstoler = 1e-4;  end
if isfield(opt,'iterbnd'), iterbnd = opt.iterbnd;  else, iterbnd = 5e1;  end

%diagonal preconditioning is used always
if isfield(opt,'verbose'), verbose = opt.verbose;  else, verbose = 1;  end


if verbose >=1
    fprintf(2,'\n<starting GNQKD> \n')
end

sizerho = size(VKVlist{1},2);
n = sizerho;
n2 = n^2;
m = length(gamma);   % number of constraints
Out = [];
bestub = log(length(VZKVlist)/length(VKVlist));
bestlb = 0;
ubs = bestub;
lbs = bestlb;
objvals = inf;
feasobjvals = inf;


const = sum_congruence(eye(size(VKVlist{1},1)),VKVlist,1) -...
    sum_congruence(eye(size(VZKVlist{1},1)),VZKVlist,1);  % constant term in the objective


MGammaorig = MGamma;
gammaorig = gamma;
pred = sqrt(diag(MGamma*MGamma'));  % left preconditioning
MGamma = diag(pred)\MGamma;
gamma = diag(pred)\gamma;

% Form a sparse null-space representation for Gamma
rp = 1:m;
[~,cp] = licols(MGamma);
ccp = setdiff(1:n2,cp);   % complement indices
cp = [cp ccp];   % permutation of cols
[~,icp] = sort(cp);   % inverse permutation of cols
MGamma = MGamma(rp,cp); % order of rows/constraints does not change problem
gamma = gamma(rp);  % for RHS
BMG = sparse(MGamma(1:m,1:m)); EMG = sparse(MGamma(1:m,m+1:end));
NGV = [sparse(BMG\EMG); -speye(n2-m)];     

MGamma_noperm = MGamma(:,icp);  % MGammas before column permutations, but with row permutations
NGV_noperm = NGV(icp,:); % matrix holding null basis for MGamma_noperm



if verbose == 2
    fprintf('\n\nin GNQKD: tolerGN %g crosstoler %g; dims:  rho %i  m %i  \n',...
        tolerGN,crosstoler,sizerho,m)
    fprintf('size MGamma %i %i and cond of MGamma*MGamma''  %g\n',...
        size(MGamma),cond(MGamma*MGamma'))
    fprintf('sparsity and cond# basis BMG  %g %g sparsity NGV %g \n',...
        nnz(BMG)/numel(BMG),condest(BMG),nnz(NGV)/numel(NGV))
end


sigmaa = 1;     % centering parameter - changes adaptively
step = 1;       % initial stepsize - changes adaptively



% start at one for initial ONE step

%%%Initializations
In = eye(n);

%Use primal feasibilities to find initial pos. def. primal vrbles
if isfield(opt,'rho0')
    rho = opt.rho0;
    rhohat = ...  % project onto feas. lin manifold
        rho - HMat(MGamma_noperm\(MGamma_noperm*Hvec(rho)-gamma));
else
    temp = MGamma\gamma;  
    rhohat = HMat(temp(icp));
    rho = In;  % Mike 
end  
d = min([eig(rhohat);0.1]);  
if d < .1   % ensure pos def
    rho = rhohat + (abs(d)+.1)*In;  % ensure pos def start
end
if min(eig(rho))<0
    rho = In;
end
pfresidual = MGammaorig*(Hvec(rhohat))-gammaorig;
if verbose == 2
    fprintf('min eig, cond #, of rho; residual of rhohat: %g %g %g \n',...
        min(eig(rho)),cond(rho),norm(pfresidual))
end
if (min(eig(rho)) > 0) && (norm(MGamma*Hvec(rho)-gamma)<tolerGN) % if initial feas point is pos. def.
    %calculate the current objective
    MDelta = sum_congruence(rho,VKVlist,0);
    MSigma = sum_congruence(rho,VZKVlist,0);
    EigDelta = eig(MDelta);
    EigSigma = eig(MSigma);
    objval = real(sum(EigDelta.*log(EigDelta)) ...
        - sum(EigSigma.*log(EigSigma)));
    objvals = [objvals objval];
    bestub = objval;
end
if isfield(opt,'y0')
    y = opt.y0;
else
    y = zeros(size(MGamma,1),1);
end
%%%%%%%%%%%%MGamma is not used again below
clear MGamma

if isfield(opt,'Z0')
    Z = opt.Z0;
else
    %evaluate gradient
    gf = sum_congruence(logm(sum_congruence_nosym(rho,VKVlist,0)),VKVlist,1)...
        +const      - ...
        sum_congruence(logm(sum_congruence_nosym(rho,VZKVlist,0)),VZKVlist,1);
    gf = (gf+gf')/2;
    d = min([eig(gf);0.1]);
    if d < .1   % ensure pos def
        Z = gf + (abs(d)+.1)*In;
        Z = (Z+Z')/2;
    else
        Z = In;
    end
end
if verbose == 2
    fprintf('min eig and cond # of initial Z: %g %g \n',...
        min(eig(Z)),cond(Z))
end
%%%%%%%%%%%%%ensure we start hermitian
gapp = real(trace(rho*Z));
mingapp = gapp;
muu = gapp/(n*2);



%calculate the primal residuals
%% initial full residual as a cell and sum up norms of els
if ~exist('gf','var')
    gf = sum_congruence(logm(sum_congruence_nosym(rho,VKVlist,0)),VKVlist,1)...
        +const      - ...
        sum_congruence(logm(sum_congruence_nosym(rho,VZKVlist,0)),VZKVlist,1);
    gf = (gf+gf')/2;
end
v = zeros(n2-m,1);   %   initial v
Fmu = {gf+HMat(MGamma_noperm'*y)-Z,...
    HMat(NGV_noperm*v)+ rhohat - rho, ... 
    Z*rho-muu*eye(n)};
Grho = sum_congruence(rho,VKVlist,0);
ZGrho = sum_congruence(rho,VZKVlist,0);
noiter = 0;
%%%%% reducing steplength if dual feas not %improv. after step=1 used
oldFmuo = norm(Fmu{1},'fro');
currFmuo = oldFmuo;
%%%%%%%%%%%%%%%%%%%%%prepare next 6 lines for evalHess
V = speye(n^2);
W = zeros(n^2);
for ii = 1:n^2
    E = HMat(V(:,ii));
    W(:,ii) = E(:);
end
W = sparse(W); 
Hess = evalHess(rho,Grho,ZGrho,VKVlist,VZKVlist,noiter,W); % hessian evaluation

Fmupr2 =  Fmu{3} + Z*Fmu{2} + ( Fmu{1} + HMat( Hess*Hvec(Fmu{2})))*rho ;  % projected residual
Fmuprnorm = norm(Fmupr2,'fro');


%calculate the current objective
MDelta = sum_congruence(rho,VKVlist,0);
MSigma = sum_congruence(rho,VZKVlist,0);
EigDelta = eig(MDelta);
EigSigma = eig(MSigma);
objval = real(sum(EigDelta.*log(EigDelta))- sum(EigSigma.*log(EigSigma)));


% Display information
VarNames = {'-log10(mingapp)','step','sigmaa','objval','pri-feas','dualfeas','cond.Jac/d','resnormaleq'};
if verbose == 2
    fprintf('%-9s','iter#')
    fprintf('%-15s',VarNames{:})
    fprintf('\n')
end


midbestublb =1 + min(mean(norm(rho,'fro')+norm(Z,'fro')), ...
                  (abs(bestub)+abs(bestlb))/2);
flagstepone = 0;  % number of times step length one taken
nodualfeasimpr = false;  % start with dual improvement happened
bestlowrho = rho;
bestlowy = y;
bestlowZ = Z;
bestlowgf = gf;
condJac = 0;  % if verbose == false
stalling = 0;   % stalling to quit while loop

%%%%%%%%%%%%start of main while loop%%%%%%%%%%%%%%%%%%%%%%%%%%%
while noiter<iterbnd && ...
     (max(mingapp,Fmuprnorm))/midbestublb > tolerGN ...
	  && (stalling <=  6 ) 
    
    noiter = noiter+1;
    ubs = [ubs bestub];
    lbs = [lbs bestlb];
    
    %%%%ensure pos. def. is preserved
    mineigs = min([eig(Grho);eig(ZGrho);eig(rho)]);
    if mineigs <= -2*eps  % machine epsilon
        stalling = stalling + 1;  % end if > 3
        fprintf('min of all eigs for logm is %g \n',mineigs)
    else
        stalling = 0; % restart to 0
    end
    
    
    %%%% Find search direction
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%change Fmu cell to vector form for lsqr/lsmr
    Fmuprvec2 = Cvec(Fmupr2);   % RHS for search direction
    % evaluate Jacobian
    Jac = zeros(2*n2,n2);
    for ii = 1:n2-m   % form Jacobian for Delta v
        NGVt = NGV_noperm(:,ii);   % extract ii-th column
        Jac(:,ii) = Cvec(Z*HMat(NGVt) + ...
            HMat(Hess*NGVt)*rho);
    end
    for ii = 1:m     % form Jacobian for Delta y
        Jac(:,n2-m+ii) = Cvec(HMat(MGamma_noperm(ii,:)')*rho);
    end
    JsJac = Jac'*Jac;
    d = sqrt(diag(JsJac));   % for diagonal preconditioning
    if verbose == 2 
        condJac = condest(diag(d)\JsJac/diag(d));
    end
    dsJac = -((Jac/diag(d))\Fmuprvec2)./d;  % diag. preconditioned soln
    deltas = dsJac;   % search direction 
    resnormeq = norm(JsJac*dsJac+Jac'*Fmuprvec2); % normal eqn residual
    
    %calculate the directions
    dv = deltas(1:end-m);
    dy = deltas(end-m+1:end);
    
    % Recover/backtrack for  drho,dZ
    if flagstepone  >= 2  % step length one already taken
        drho = HMat(NGV_noperm*dv);
    else
        drho = HMat(NGV_noperm*dv)+Fmu{2};
    end
    dZ =  HMat(Hess*Hvec(drho))+HMat(MGamma_noperm'*dy) + Fmu{1};
    %%%% search direction found
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% Find step length; check the semidefiniteness of solution
    rhop = rho + step*drho;
    [~,choltruerho] = chol(rhop);
    while choltruerho  > 0 && step > 1e-9
        step = .8*step;  % original backtracking
        rhop = rho + step*drho;
        [~,choltruerho] = chol(rhop);
    end
    Zp = Z + step*dZ;
    [~,choltrueZ] = chol(Zp);
    while  choltrueZ > 0 && step > 1e-9
        step = .8*step;  % original backtracking
        Zp = Z + step*dZ;
        [~,choltrueZ] = chol(Zp);
    end
    rhop = rho + step*drho;
    %%%%%%%%%%%% add backtracking for obj function
    [~,choltrueGrho] = chol(sum_congruence(rhop,VKVlist,0));
    while choltrueGrho  > 0 && step > 1e-9
        step = .8*step;  % original backtracking
        rhop = rho + step*drho;
        [~,choltrueGrho] = chol(sum_congruence(rhop,VKVlist,0));
    end
    %%%% following is debug check 
    if min(eig(sum_congruence(rhop,VKVlist,0)))<= -2*eps
        fprintf('min eig Grho <= -2*eps \n')
        keyboard
    end
    [~,choltrueZGrho] = chol(sum_congruence(rhop,VZKVlist,0));
    while choltrueZGrho  > 0 && step > 1e-9
        step = .8*step;  % original backtracking
        rhop = rho + step*drho;
        [~,choltrueZGrho] = chol(sum_congruence(rhop,VZKVlist,0));
    end
    %%%% following is debug check
    if min(eig(sum_congruence(rhop,VZKVlist,0)))<= -2*eps
        fprintf('min eig ZGrho <= -2*eps \n')
        keyboard
    end
    
    if step > (1/.97) && flagstepone == 0  % Newton step = 1 not taken yet
        if verbose == 2 
            fprintf(2,'FIRST Newton step = 1 taken\n');
        end
        flagstepone = 1;  % indicates that step = 1 taken for future steps
        step = 1;  % take a Newton step=1  exactly since first time
    elseif step > 1e-9
        %The following is done ONLY the iteration after the first step=1
        if flagstepone == 1 % Newton step taken last iteration
            if verbose == 2
                fprintf(2,'exact pf starts; setting flagstepone = 2 \n');
            end
            flagstepone = 2;
        end
        step = .97*step;   % backtrack
    else
	    fprintf('In GNQKD, STOP: step lengths are too small \n')
	    break  % get out of main while loop
    end
    %%%% heuristic for dualfeasibility
    if flagstepone >= 2  && nodualfeasimpr == true
        step = .5*step;  % backtrack if no dualfeas impr.
    end
    %%%% steplength found; update variables
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
    
    %%%% variables/residuals update
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    v = v + step*dv;  % update v for simple/exact rho update
    % by linearity  newrho = rhohat+N^*v+step*N^*dv
    if flagstepone >= 2   % step = 1 taken at least one iter previously
        rho = rhohat + HMat(NGV_noperm*v); % exact primal feas.
        Fmu{2} = zeros(n);

    else
        rho = rho + step*drho;
        Fmu{2} = HMat(NGV_noperm*v)+ rhohat - rho;
    end
    y = y + step*dy;
    Z = Z + step*dZ;

    Grho = sum_congruence(rho,VKVlist,0);
    ZGrho = sum_congruence(rho,VZKVlist,0);
    
    mineigs = min([eig(Grho);eig(ZGrho);eig(rho)]);
    if mineigs <= -2*eps  % machine epsilon
        stalling = stalling + 1;  % end if > 3
        fprintf('min of all eigs for logm is %g \n',mineigs)
    else
        stalling = 0; % restart to 0
    end
    gf = sum_congruence(logm(sum_congruence_nosym(rho,VKVlist,0)),VKVlist,1)...
        +const      - ...
        sum_congruence(logm(sum_congruence_nosym(rho,VZKVlist,0)),VZKVlist,1);
    Fmu{1} = gf+HMat(MGamma_noperm'*y)-Z;
    temp = norm(Fmu{1},'fro');
    if temp > 0   % update last two positive dual residuals
        oldFmuo = currFmuo;
        currFmuo = temp;
        if currFmuo > oldFmuo
            nodualfeasimpr = true;  % no improvement in dualfeas
        else
            nodualfeasimpr = false;
        end
    end
    Fmu{3} = Z*rho-muu*eye(n);
    
    % Hessain evaluated after taking a new step
    Hess = evalHess(rho,Grho,ZGrho,VKVlist,VZKVlist,noiter,W);   
    
    if flagstepone >= 2  % exact primal feas. obtained
        Fmupr2 = Fmu{3} + Z*Fmu{2} +  Fmu{1}*rho ;
    else
        Fmupr2 =  ...
            Fmu{3} + Z*Fmu{2} + ( Fmu{1} + HMat( Hess*Hvec(Fmu{2})))*rho ;
    end
    
    %%%% End of variables/residuals update
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %calculate the current objective
    MDelta = sum_congruence(rho,VKVlist,0);
    MSigma = sum_congruence(rho,VZKVlist,0);
    EigDelta = eig(MDelta);
    EigSigma = eig(MSigma);
    objval = real(sum(EigDelta.*log(EigDelta))- sum(EigSigma.*log(EigSigma)));
    
    Fmuprnorm = norm(Fmupr2,'fro');
    gapp = real(trace(rho*Z));
    mingapp = min(gapp,bestub-bestlb);
    midbestublb =1 + min(mean(norm(rho,'fro')+norm(Z,'fro')), ...
                  (abs(bestub)+abs(bestlb))/2);
              
    %% if sufficient accuracy test for valid upper/lower bnds
    if       (Fmuprnorm/(abs(objval)+1) < crosstoler) || ...
		    mod(noiter,10) == 0
        %%%%%%%%%%%check for ub %%%%%%%%%%%%%%%%%%%
        if flagstepone < 1   % step length one NOT taken yet
            newrho = ...  % project onto feas. lin manifold
                rho - HMat(MGamma_noperm\(MGamma_noperm*Hvec(rho)-gamma));
        else
            newrho = rho;
        end
        if min(eig(newrho))>0  % test psd feas.
            ub = ...   % valid upper bnd
                real(trace(sum_congruence(newrho,VKVlist,0)*logm(sum_congruence_nosym(newrho,VKVlist,0))))-...
                real(trace(sum_congruence(newrho,VZKVlist,0)*logm(sum_congruence_nosym(newrho,VZKVlist,0))));
            ubs = [ubs ub];
            bestub = min(bestub,ub);
        end
        %%%%%%%%%%%end check for ub %%%%%%%%%%%%%%%%%%%
        
        %%%%%%%%%%%check for lower bound 
        dualbdZ =  HMat(MGamma_noperm'*y);
	%%%%check Farkas lemma alternative for infeas of primal
	if gamma'*y < -1e-12  && min(eig(dualbdZ)) >=0 
		fprintf('inGNQKD: primal infeas detected\n')
		keyboard
		break
	end
	%% continue with check on lower bound
        newZ =  gf+dualbdZ;  % nearest dual feas.
        mineignewZ = min(eig(newZ));
        if mineignewZ >0   % if psd feas. get valid lower bnd
            if flagstepone >= 2  % step one taken already
                lb = real(objval - trace(rho*newZ));
            else
                lb = real(objval + y'*(MGamma_noperm*Hvec(rho)-gamma) ...
                    - real(trace(rho*newZ)));
            end
            lbs = [lbs lb];
            if lb >= bestlb
                bestlb = max(bestlb,lb);
                bestlowrho = rho;
                bestlowy = y;
                bestlowZ = Z;
                bestlowgf = gf;
            end
            %%%% test for 'little' loss in compl slack
            if real(trace(newZ*rho)) <= 2*gapp && ...
                    mineignewZ >= min(eig(Z))
                Z = newZ;   % use improved dual feas.
                if verbose == 2 
                    fprintf(2,'LB updated:  bestub/bestlb/rel-gap: ')
                    fprintf(2,'%g %g %g ',...
                        bestub,bestlb,...
                        (bestub-bestlb)/midbestublb);
                    fprintf('AND: newZ used for dual feasibility\n')
                end
                %%% update residual changes
                Fmu{1} = zeros(n);  
                Fmupr2 = Fmu{3} + Z*Fmu{2} + ...
                         ( HMat( Hess*Hvec(Fmu{2})))*rho ;
            else
                if verbose == 2
                    fprintf(2,'LB updated:  bestub/bestlb/rel-gap: %g %g %g \n',...
                        bestub,bestlb,...
                        (bestub-bestlb)...
                        /midbestublb);
                end
            end
        end
    end
    muu = sigmaa*gapp/(2*n);
    
    if verbose == 2 
        fprintf('%-9i',noiter)
        fprintf('%-15e',-log10((mingapp)/(abs(objval)+1)),step,sigmaa,...
            objval,norm(Fmu{2}),norm(Fmu{1},'fro'),condJac,resnormeq);
        fprintf('\n')
    end
    gapps(noiter) = gapp;
    mingapps(noiter) = mingapp;
    steps(noiter) = step;
    sigmaas(noiter) = sigmaa;
    objvals(noiter) = objval;
    condJacs(noiter) = condJac;
    resnormeqs(noiter) = resnormeq;
    
    
    % update muu and sigmaa (centering parameter) depending on stepsize
    if step < .001
        muu = muu * 2.1;
        sigmaa = 1;
        step = .9;  %muu increased
    elseif step < .01
        muu = muu * 2;
        sigmaa = 1;
        step = 1.01;  % muu unchanged
    elseif step < .5
        muu = muu * 1.5;
        sigmaa = .5;
        step = 1.01;
    elseif step < .95
        muu = muu * 1.2;
        sigmaa = (1-.3*min(step,1));
        step = 1.2;
    elseif step < 1.05
        muu = muu *.5;
        sigmaa = min((1-.3*min(step,1))/2,...
            norm(Z,'fro'));
        step = 1.5;
    else
        sigmaa = min(1e-3,norm(Z,'fro'));
        step = 1.5;
    end
end          % end of main while loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose == 2 
    fprintf('after WHILE: noiter; LHStolerGN; stalling  %i  %g  %i \n',...
	    noiter,(max(mingapp,Fmuprnorm))/midbestublb,stalling)
end

    if length(lbs) < noiter
        lbs(length(lbs)+1:noiter) = lbs(end);
    end
    if length(ubs) < noiter
        ubs(length(ubs)+1:noiter) = ubs(end);
    end
    if verbose >= 1
        fprintf('%-9s','iter#')
        fprintf('%-15s',VarNames{:})
        fprintf('\n')
        fprintf('%-9i',noiter)
        fprintf('%-15e',-log10((mingapp)/(abs(objval)+1)),step,sigmaa,...
            objval,norm(Fmu{2}),norm(Fmu{1},'fro'),condJac,resnormeq);
        fprintf('\n')
        if noiter < iterbnd && stalling <= 6
            fprintf('USER: from GNQKD.m: Toler achieved - algorithm ends\n')
        end
        if noiter >= iterbnd  && stalling <= 6
            fprintf('USER: from GNQKD.m: ERROR;  exceeded iter bnd: %i  \n',...
                noiter);
        end
        if stalling > 6
            fprintf('USER: from GNQKD.m: ERROR;  stopping due to stalling')
        end
    end
    
    
    if verbose == 2 
        fprintf('[mineig, condest]: rho,Z,delta,sigma, ')
        fprintf('[%g %g]; [%g %g];  [%g %g]; [%g %g] \n',...
            min(eig(rho)),condest(rho),...
            min(eig(Z)),condest(Z),...
            min(eig(MDelta)),condest(MDelta),...
            min(eig(MSigma)),condest(MSigma))
    end
    
    if verbose >= 1 
        fprintf(2,'FINAL in GNQKD:  ub/lb/rel-gap/rel||RES||: ')
        fprintf('%g %g %g %g \n', ...
                 bestub,bestlb, ...
                (bestub-bestlb)/midbestublb, ...
                 Fmuprnorm/midbestublb     );
    end
    Out.mingapps = mingapps;
    Out.rho = rho;
    Out.y = y;
    Out.Z = Z;
    Out.muu = muu;
    Out.objval = bestub;
    Out.objval2 = bestub/log(2);
    Out.bestub = bestub;
    Out.bestlb = bestlb;
    Out.ubs = ubs;
    Out.lbs = lbs;
    Out.noiter = noiter;
    Out.gapps = gapps;
    Out.steps = steps;
    Out.sigmaas = sigmaas;
    Out.objvals = objvals;
    Out.condJacs = condJacs;
    Out.resnormeqs = resnormeqs;
    Out.bestlowrho = bestlowrho;  
    Out.bestlowy = bestlowy;      
    Out.bestlowZ = bestlowZ;     
    Out.bestlowgf = bestlowgf;    
    Out.gradf = gf;
    if strcmp(profilef,'solver')
        fprintf('Profile at end of solver file\n')
        profile report
        keyboard
    end
    
    if verbose >= 1
        fprintf(2,'<end of GNQKD> \n')
    end
end   % end of function

