sedumi.m

经济学专业代码

    r0 = sum(R.w);
    cx = by + y0*R.sd - x0 / d.l(1);
    rgap = max(cx-by,0) / max([abs(cx),abs(by),1e-3 * x0]);
    precision1=y0*r0/(1+x0);
    precision2=(y0 * r0 + rgap)/x0;
    my_fprintf(pars.fid,'%1.1E\n',max(precision1,precision2));
    if precision1 < pars.eps       % P/D residuals small
        if precision2 < pars.eps    %Approx feasible and optimal
            STOP = 1;
            break
        elseif y0 * R.maxRb + x0 * R.maxb < -pars.eps * cx   % Approx Farkas
            STOP = 1;
            break
        elseif y0 * R.maxRc + x0 * R.maxc < pars.eps * by    % Approx Farkas
            STOP = 1;
            break;
        end
    end
    if iter >= pars.maxiter
        my_fprintf(pars.fid,'Maximum number of iterations reached.\n');
        STOP = -1;
    end
end % while STOP == 0.
my_fprintf(pars.fid,'\n');
clear ADA 
nnzLadd=nnz(L.add);
nnzLskip=nnz(L.skip);
normLL=full(max(max(abs(L.L))));
clear L
% ************************************************************
% FINAL TASKS:
% ************************************************************
cputime2=cputime;
info.iter = iter;
info.feasratio = feasratio;
info.pinf = 0; info.dinf = 0;
info.numerr = 0;
% ------------------------------------------------------------
% Create x = D*v.
% ------------------------------------------------------------
if STOP == 2                % Exact optimal solution found (LP)
    x = xsol;
    y = ysol;
elseif STOP == 3            % Farkas solution y found (in early stage)
    x = zeros(length(c),1);
else
    x = [sqrt(d.l).*v(1:K.l); asmDxq(d,v,K); psdscale(d,v,K,1)];
end
% --------------------------------------------------
% Compute cx, Ax, etc.
% --------------------------------------------------
x0 = x(1);
cx = full(sum(c.*x)); 
abscx = sum(abs(c).*abs(x));
by = full(sum(b.*y));
Ax = Amul(A,dense,x,0);
Ay = full(Amul(A,dense,y,1));      % "full" since y may be scalar.
normy = norm(y);
normx = norm(x(2:end));
clear A
% ------------------------------------------------------------
% Determine infeasibility
% ------------------------------------------------------------
pinf = norm(x0*b-Ax);
z = qreshape(Ay-x0*c,1,K);
dinf = max(eigK(z,K));
if x0 > 0
    relinf = max(pinf / (1+R.maxb), dinf / (1+R.maxc)) / x0;
    % ------------------------------------------------------------
    % If infeasibility larger than epsilon, evaluate Farkas-infeasibility
    % ------------------------------------------------------------
    if relinf > pars.eps
        pdirinf = norm(Ax);
        ddirinf = max(eigK(qreshape(Ay,1,K),K));
        if cx < 0.0
            reldirinf = pdirinf / (-cx);
        else
            reldirinf = inf;
        end
        if by > 0.0
            reldirinf = min(reldirinf, ddirinf / by);
        end
        % ------------------------------------------------------------
        % If the quality of the Farkas solution is good and better than
        % the approx. feasible soln, set x0=0: Farkas solution found.
        % ------------------------------------------------------------
        if (reldirinf < pars.eps) | (relinf > max(pars.bigeps, reldirinf))
            x0 = 0.0;
            pinf = pdirinf;
            dinf = ddirinf;
        end
    end % relinf too large
end % x0 > 0
% ------------------------------------------------------------
% Interpret the solution as feasible:
% ------------------------------------------------------------
if x0 > 0
    x = x / x0;
    y = y / x0;
    pinf = pinf /x0;
    dinf = dinf / x0;
    cx = cx/ x0;
    by = by / x0;
    normx = normx / x0;
    normy = normy / x0;
    if cx <= by                % zero or negative duality gap
        sigdig = Inf;
    elseif cx == 0.0           % Dual feasibility problem
        sigdig = -log(-by/(R.maxb*normy +1E-10 * x0)) / log(10);
    elseif by == 0.0           % Primal feasibility problem
        sigdig = -log(cx / (R.maxc*normx +1E-10 * x0)) / log(10);
    else                       % Optimization problem
        sigdig = -log((cx-by)/(abs(by) + 1E-5 * (x0+abscx))) / log(10);
    end
    my_fprintf(pars.fid,...
        'iter seconds digits       c*x               b*y\n');
    my_fprintf(pars.fid,'%3d %8.1f %5.1f %- 17.10e %- 17.10e\n',...
        iter,cputime2-cputime1,sigdig,cx,by);
    my_fprintf(pars.fid,'|Ax-b| = %9.1e, [Ay-c]_+ = %9.1E, |x|=%9.1e, |y|=%9.1e\n',...
        pinf,dinf,normx,normy);
    % ------------------------------------------------------------
    % Determine level of numerical problems with x0>0 (feasible)
    % ------------------------------------------------------------
    if STOP == -1
        r0 = max([10^(-sigdig); pinf;dinf] ./ [1; 1+R.maxb+(1E-3)*R.maxRb;...
            1+R.maxc+(1E-3)*R.maxRc]);
        if r0 > pars.bigeps
            my_fprintf(pars.fid, 'No sensible solution found.\n');
            info.numerr = 2;                          % serious numerical error
        elseif r0 > pars.eps
            info.numerr = 1;                          % moderate numerical error
        else
            info.numerr = 0;                          % achieved desired accuracy
        end
    end
else  % (if x0>0)
    % --------------------------------------------------
    % Infeasible problems: pinf==norm(Ax), dinf==max(eigK(At*y,K)).
    % --------------------------------------------------
    if pinf < -pars.bigeps * cx
        info.dinf = 1;
        abscx = -cx;
        pinf = pinf / abscx;
        normx = normx / abscx;
        x = x / abscx;
        my_fprintf(pars.fid, 'Dual infeasible, primal improving direction found.\n');
    end
    if dinf < pars.bigeps * by
        info.pinf = 1;
        dinf = dinf / by;
        normy = normy / by;
        y = y / by;
        my_fprintf(pars.fid, 'Primal infeasible, dual improving direction found.\n');
    end
    my_fprintf(pars.fid,'iter seconds  |Ax|    [Ay]_+     |x|       |y|\n');
    my_fprintf(pars.fid,'%3d %8.1f %9.1e %9.1e %9.1e %9.1e\n',...
        iter,cputime2-cputime1,pinf,dinf,normx,normy);
    % --------------------------------------------------
    % Guess infeasible, but stopped due to numerical problems
    % --------------------------------------------------
    if info.pinf + info.dinf == 0
        my_fprintf(pars.fid, 'Failed: no sensible solution/direction found.\n');
        info.numerr = 2;
    elseif STOP == -1
        if (pinf > -pars.eps * cx) & (dinf > pars.eps * by)
            info.numerr = 1;
        else
            info.numerr = 0;
        end
    end
end
% ----------------------------------------
% - Bring xsol into the complex format of original (At,c),
% - transform q-variables into r-variables (Lorentz),
% - bring ysol into complex format, indicated by K.ycomplex.
% - at 0's in ysol where rows where removed.
% ----------------------------------------'
[x,y,K] = posttransfo(x,y,prep,K,pars);
% Detailed timing
%Preprocessing+IPM+Postprocessing
info.timing=[cputime1-cputime0 cputime2-cputime1 cputime-cputime2];
% Total time (for backward compatibility)
info.cpusec=sum(info.timing);
my_fprintf(pars.fid,'\nDetailed timing (sec)\n')
my_fprintf(pars.fid,'   Pre          IPM          Post\n')
my_fprintf(pars.fid,'%1.3E    ',info.timing)
my_fprintf(pars.fid,'\n')


% ----------------------------------------
% Make a fancy v-plot if desired
% ----------------------------------------
if pars.vplot == 1
    subplot(2,1,1)
    plot(0:iter,vlist,'o',[0 iter],[1 1],'b',...
        [0 iter],[pars.theta pars.theta],'g')
    title('Wide region v-plot')
    xlabel('iterations')
    ylabel('normalized v-values')
    subplot(2,1,2)
    plot(1:iter,ratelist)
    axis([0 iter 0 1])
    title('Reduction rates')
    xlabel('iterations')
    ylabel('reduction rate')
end
my_fprintf(pars.fid,'Max-norms: ||b||=%d, ||c|| = %d,\n',R.maxb,R.maxc);
my_fprintf(pars.fid,'Cholesky |add|=%d, |skip| = %d, ||L.L|| = %g.\n',...
    nnzLadd, nnzLskip, normLL);

% ----------------------------------------
% Compute error measures if needed
% ----------------------------------------
if ~isempty(origcoeff)
    %Reload the original coefficients
    s=(origcoeff.c)-(origcoeff.At)*sparse(y);       %To make s sparse
    cx=sum((origcoeff.c).*x);        %faster than c'*x
    by=sum((origcoeff.b).*y);
    xs=sum(x.*s);
    normb=norm(origcoeff.b,1);
    normc=norm(origcoeff.c,1);
    info.err=zeros(1,6);
    %     Error measures.
    %     Primal infeasibility
    info.err(1)=norm(x'*(origcoeff.At)-(origcoeff.b)',2)/(1+normb);
    %Let us get rid of the K.f part, since the free variables don't make
    %any difference in the cone infeasibility.
    %origcoeff.K.f=0;
    %     Primal cone infeasibility
    info.err(2)=max(0,-min(eigK(full(x(origcoeff.K.f+1:end)),origcoeff.K))/(1+normb));
    %     Dual infeasibility
    %info.err(3)=0.0; %s is not maintained explicitely
    %     Dual cone infeasibility
    info.err(4)=max(0,-min(eigK(full(s(origcoeff.K.f+1:end)),origcoeff.K))/(1+normc));
    %     Relative duality gap
    info.err(5)=(cx-by)/(1+abs(cx)+abs(by));
    %     Relative complementarity
    info.err(6)=xs/(1+abs(cx)+abs(by));

    my_fprintf(pars.fid,'\nDIMACS error measures\n')
    my_fprintf(pars.fid,'  PInf     PConInf     DInf    DConInf    RelGap    RelComp\n')
    my_fprintf(pars.fid,'%2.2E  ',info.err)
    my_fprintf(pars.fid,'\n')
end