sedumi.m

经济学专业代码

    end
end
% ----------------------------------------
% Bring data in (strict) internal SeDuMi form, check parameters,
% and print welcome.
% ----------------------------------------
[A,b,c,K,prep,origcoeff] = pretransfo(A,b,c,K,pars);
if prep.cpx.dim>0
    origcoeff=[];      % No error measures for complex problems.
end
lponly = (K.l == length(c));
pars = checkpars(pars,lponly);
% ----------------------------------------
% Print welcome and statistics of cone-problem
% ----------------------------------------
my_fprintf(pars.fid,'SeDuMi 1.1R3 by AdvOL, 2006 and Jos F. Sturm, 1998-2003.\n');
switch pars.alg
    case 0
        my_fprintf(pars.fid,'Alg = 0: No corrector, ');
    case  1
        my_fprintf(pars.fid,'Alg = 1: v-corrector, ');
    case 2
        my_fprintf(pars.fid,'Alg = 2: xz-corrector, ');
end
switch pars.stepdif
    case 0
    case 1
        my_fprintf(pars.fid,'Step-Differentiation, ');
    case 2
        my_fprintf(pars.fid,'Adaptive Step-Differentiation, ');
end
my_fprintf(pars.fid,'theta = %5.3f, beta = %5.3f\n',pars.theta,pars.beta);
% --------------------------------------------------
% Print preprocessing information
% --------------------------------------------------
if pars.prep==1
    if isfield(prep,'sdp')
        blockcount=0;
        varcount=0;
        for sdpind=1:length(prep.sdp)
            if prep.sdp{sdpind}(1)==1
                blockcount=blockcount+1;
                varcount=varcount+prep.sdp{sdpind}(2);
            end
        end
        if blockcount>0
            my_fprintf(pars.fid,'Detected %i diagonal SDP block(s) with %i linear variables\n',blockcount,varcount);
        end
    end
    if isfield(prep,'freeblock1') & length(prep.freeblock1)>0
        my_fprintf(pars.fid,'Detected %i free variables in the linear part\n',length(prep.freeblock1));
    end
    if isfield(prep,'Kf') & prep.Kf>0
        switch pars.free
            case 0
                my_fprintf(pars.fid,'Split %i free variables\n',prep.Kf);
            case 1
                my_fprintf(pars.fid,'Put %i free variables in a quadratic cone\n',prep.Kf);
        end
    end
end
% --------------------------------------------------
% Remove dense columns (if any)
% --------------------------------------------------
Ablkjc = partitA(A,K.mainblks);
[dense,DAt.denq] = getdense(A,Ablkjc,K,pars);
if length(dense.cols) > 0
    dense.A = A(dense.cols,:)';
    A(dense.cols,:) = 0.0;
    Ablkjc = partitA(A,K.mainblks);
else
    dense.A = sparse(length(b),0);
end
% ----------------------------------------
% Order constraints from sparse to dense, and find corresponding
% incremental nonzero pattern "Aord.dz" of At*dy in PSD part.
% ----------------------------------------
Aord.lqperm = sortnnz(A,[],Ablkjc(:,3));        % Sparse LP+Lorentz
DAt.q = findblks(A,Ablkjc,2,3,K.qblkstart);     % Lorentz ddotA-part
if ~isempty(DAt.q)
    DAt.q(dense.q,:) = 0.0;
    DAt.q = DAt.q + spones(extractA(A,Ablkjc,1,2,K.mainblks(1),K.mainblks(2)));
    Aord.qperm = sortnnz(DAt.q,[],[]);
else
    Aord.qperm = (1:length(b))';
end
[Aord.sperm, Aord.dz] = incorder(A,Ablkjc(:,3),K.mainblks(3));  % PSD
% ----------------------------------------
% Get nz-pattern of ADA.
% ----------------------------------------
ADA = getsymbada(A,Ablkjc,DAt,K.sblkstart);
% ----------------------------------------
% Ordering and symbolic factorization of ADA.
% ----------------------------------------
L = symbchol(ADA);
% --------------------------------------------------
% Symbolic fwsolve dense cols: L\[dense.A, dense.blkq],
% sparse ordering for dense column factorization
% --------------------------------------------------
symLden = symbcholden(L,dense,DAt);
% ----------------------------------------
% Initial solution
% ----------------------------------------
[d, v,vfrm, y,y0, R] = sdinit(A,b,c,dense,K,pars);
n = length(vfrm.lab);                         % order of K
merit = (sum(R.w) + max(R.sd,0))^2 * y0 / R.b0;  % Merit function
my_fprintf(pars.fid,'eqs m = %g, order n = %g, dim = %g, blocks = %g\n',...
    length(b),n,length(c),1 + length(K.q) + length(K.s));
my_fprintf(pars.fid,'nnz(A) = %d + %d, nnz(ADA) = %d, nnz(L) = %d\n',nnz(A),nnz(dense.A), nnz(ADA), nnz(L.L));
if length(dense.cols) > 0
    my_fprintf(pars.fid,'Handling %d + %d dense columns.\n',...
        length(dense.cols),length(dense.q));
end
my_fprintf(pars.fid,' it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec\n');
my_fprintf(pars.fid,'  0 :            %8.2E %5.3f\n',merit,0);
% ----------------------------------------
% Initialize iterative statistics
% ----------------------------------------
STOP = 0;
err.maxb = 0.0;                 % Estimates the need to recompute residuals
iter = 0;
if pars.vplot == 1
    vlist = [vfrm.lab];
    ratelist = [];
end
wr.delta = 0.0;
wr.desc = 1;
%Seems unnecessary
%rate = 1.0;
feasratio = 0.0;
cputime1 = cputime;
% ************************************************************
% MAIN PREDICTOR-CORRECTOR LOOP
% ************************************************************
while STOP == 0
    iter = iter+1;
    if any(iter == pars.stopat)
        keyboard
    end
    
    if pars.stepdif==2 & ...
            (iter>20 | (iter>1 & (err.kcg + Lsd.kcg>3)) | ...
            (iter>5 & abs(1-feasratio)<0.05) )
        pars.stepdif=1;
    end
    % --------------------------------------------------
    % Compute ADA
    % --------------------------------------------------
    DAt = getDAtm(A, Ablkjc, dense, DAt.denq, d, K);
    ADA = getada1(ADA, A, Ablkjc(:,3), Aord.lqperm, d, K.qblkstart);
    ADA = getada2(ADA, DAt, Aord, K);
    [ADA,absd] = getada3(ADA, A, Ablkjc(:,3), Aord, invcholfac(d.u, K, d.perm), K);
    % ------------------------------------------------------------
    % Block Sparse Cholesky: ADA(L.perm,L.perm) = L.L*diag(L.d)*L.L'
    % ------------------------------------------------------------
    [L.L,L.d,L.skip,L.add] = blkchol(L,ADA,pars.chol,absd);
    % ------------------------------------------------------------
    % Factor dense columns
    % ------------------------------------------------------------
    [Lden, L.d] = deninfac(symLden, L,dense,DAt,d,absd, K.qblkstart,pars.chol);
    % ----------------------------------------
    % FACTORIZATION of self-dual embedding
    % ----------------------------------------
    Lsd = sdfactor(L,Lden, dense,DAt, d,v,y, A,c,K,R, y0,pars);
    % ------------------------------------------------------------
    % Compute and take IPM-step
    % from (v,y,v, y0) --> (xscl,y,zscl,y0)
    % ------------------------------------------------------------
    y0Old = y0;
    [xscl,yNxt,zscl,y0Nxt, w,relt, dxmdz,err, wr] = ...
        wregion(L,Lden,Lsd,...
        d,v,vfrm,A,DAt,dense, R,K,y,y0,b, pars, wr);
    % ------------------------------------------------------------
    % Evaluate the computed step.
    % ------------------------------------------------------------
    if y0Nxt > 0
        R.b = R.b + err.b / y0Nxt;
        R.sd = R.sd + err.g / y0Nxt;
        R.b0 = R.b0 + err.db0 / y0Nxt;
        y0 = y0Nxt;
    else
        R.b = (y0Nxt * R.b + err.b) / y0Old;      % In fact, we should have y0=0.
        R.sd = (y0Nxt * R.sd + err.g) / y0Old;
        R.b0 = (y0Nxt * R.b0 + err.db0) / y0Old;
        R.w(2) = abs(y0Nxt/y0Old) * R.w(2);        %=0: dual feasible
        R.c = (y0Nxt/y0Old) * R.c;
        R.maxRc = norm(R.c,inf);
        y0 = y0Old;
    end
    R.maxRb = norm(R.b,inf);                % Primal residual
    R.w(1) = 2 * pars.w(1) * R.maxRb / (1+R.maxb);
    meritOld = merit;
    merit = (sum(R.w) + max(R.sd,0))^2 * y0 / R.b0;
    rate = merit / meritOld;
    if (rate >= 0.9999) & (wr.desc == 1)
        % ------------------------------------------------------------
        % STOP = -1  --> Stop due to numerical problems
        % ------------------------------------------------------------
        STOP = -1;                  % insuf. progress in descent direction.
        iter = iter - 1;
        y0 = y0Old;
        my_fprintf(pars.fid,'Run into numerical problems.\n');
        break
    end
    feasratio = dxmdz(1) / v(1);            % (deltax0/x0) - (deltaz0/z0)
    % --------------------------------------------------
    % Primal-Dual transformation
    % --------------------------------------------------
    y = yNxt;
    by = full(sum(b.*y));
    [d,vfrm] = updtransfo(xscl,zscl,w,d,K);
    v = frameit(vfrm.lab,vfrm.q,vfrm.s,K);
    x0 = sqrt(d.l(1)) * v(1);
    % ----------------------------------------
    % SHOW ITERATION STATISTICS
    % ----------------------------------------
    my_fprintf(pars.fid,' %2.0f : %10.2E %8.2E %5.3f %6.4f %6.4f %6.4f %6.2f %2d %2d  ',...
        iter,by/x0,merit,wr.delta,rate,relt.p,relt.d,feasratio,err.kcg, Lsd.kcg);
    if pars.vplot == 1
        vlist = [vlist vfrm.lab/sqrt((R.b0*y0)/n)];
        ratelist = [ratelist rate];
    end
    % ----------------------------------------
    % If we get in superlinear region of LP,
    % try to guess optimal solution:
    % ----------------------------------------
    if lponly & (rate < 0.05)
        [xsol,ysol] = optstep(A,b,c, y0,y,d,v,dxmdz, ...
            K,L,symLden,dense, Ablkjc,Aord,ADA,DAt, feasratio, R,pars);
        if ~isempty(xsol)
            STOP = 2;                   % Means that we guessed right !!
            feasratio = 1 - 2*(xsol(1)==0);
            break
        end
    elseif (by > 0) & (abs(1+feasratio) < 0.05) & (R.b0*y0 < 0.5)
        if max(eigK(full(qreshape(Amul(A,dense,y,1),1,K)),K)) <= pars.eps * by
            STOP = 3;                   % Means Farkas solution found !
            break
        end
    end
    % --------------------------------------------------
    % OPTIMALITY CHECK: stop if y0*resid < eps * (x0+z0).
    % For feas. probs, we should divide the residual by x0, otherwise by z0.
    % Before stopping, recompute R.norm, since it may have changed due to
    % residual updates (the change should be small though).
    % --------------------------------------------------