diff --git a/.gitignore b/.gitignore
index e05e1d1..f0fe8de 100644
--- a/.gitignore
+++ b/.gitignore
@@ -216,6 +216,13 @@ pythontex-files-*/
TSWLatexianTemp*
## Editors:
+
+# MATLAB editor
+*.asv
+
+# gedit
+*.*~
+
# WinEdt
*.bak
*.sav
@@ -228,5 +235,3 @@ TSWLatexianTemp*
# KBibTeX
*~[0-9]*
-admmdual/dataChorDec.m
-admmdual/dataChorDec.m
diff --git a/README.md b/README.md
index fca05d6..79fa15e 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,15 @@
# CDCS
CDCS (Cone Decomposition Conic Solver) is an open-source MATLAB solver for sparse conic programs with partially decomposable conic constraints. CDCS implements the alternating direction method of multipliers (ADMM)
-described in on our paper [_Fast ADMM for Semidefinite Programs with Chordal Sparsity_](https://arxiv.org/pdf/1609.06068v2.pdf) (included in the `doc/` folder).
+described in our papers [_Fast ADMM for Semidefinite Programs with Chordal Sparsity_](https://arxiv.org/pdf/1609.06068v2.pdf) and [_Fast ADMM for homogeneous self-dual embeddings of sparse SDPs_](https://arxiv.org/pdf/1611.01828.pdf) (included in the `doc/` folder)
-**Current version:** 1.0
+**Current version:** 1.1
-**Release notes:** CDCS is based on a temporary research code called ADMM-PDCP, which is no longer maintained. If you downloaded ADMM-PDCP, please replace it with CDCS.
+**Release notes:**
+
+* Homogeneous self-dual embedding is the new default method
+* The termination codes have changed. This means that if you use CDCS from YALMIP, the termination code returned by YALMIP will be incorrect. This should be fixed in the next YALMIP release!
+* CDCS is based on a temporary research code called ADMM-PDCP, which is no longer maintained. If you downloaded ADMM-PDCP, please replace it with CDCS.
## Contents
@@ -21,7 +25,7 @@ described in on our paper [_Fast ADMM for Semidefinite Programs with Chordal Spa
CDCS solves in the standard primal and dual vectorized forms
minimize c'x maximize b'y
- (1) subject to Ax = b, (2) subject to c - A'y = z,
+ (1) subject to Ax = b, (2) subject to A'y + z = c,
x \in K z \in K*
where the conic constraint `x \in K` are partially decomposable. This means that
@@ -41,6 +45,9 @@ sparsity pattern. The other supported cone types are not decomposed.
This means that CDCS is most suitable for large sparse semidefinite programs (SDPs),
although it can be used for any conic program over the supported cones.
+CDCS offers a choice to solve the primal problem (1) only, the dual problem (2) only,
+or the homogeneous self-dual embedding of the two problems. _From version 1.1, the
+homogeneous self-dual embedding is the default method._
## Quick start
@@ -72,26 +79,36 @@ you have any suggestions for improvement, or find any bugs, feel free to [contac
## How to cite
-If you find CDCS useful, please cite:
+If you find CDCS useful, please cite at least one of the following papers as appropriate:
```
-@article{ZFPGW2016,
- archivePrefix = {arXiv},
- eprint = {1609.06068v2},
- primaryClass = "math-OC",
- author = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
- title = {{Fast ADMM for Semidefinite Programs with Chordal Sparsity}}
- }
+@article{{ZFPGWhsde2016,
+ archivePrefix= {arXiv},
+ eprint = {1611.01828},
+ primaryClass = "math-OC",
+ author = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
+ title = {{Fast ADMM for homogeneous self-dual embeddings of sparse SDPs}}
+ }
+
+@article{ZFPGWpd2016,
+ archivePrefix= {arXiv},
+ eprint = {1609.06068v2},
+ primaryClass = "math-OC",
+ author = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
+ title = {{Fast ADMM for Semidefinite Programs with Chordal Sparsity}}
+ }
@misc{CDCS,
author = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
- title = {{CDCS}: Cone Decomposition Conic Solver, version 1.0},
+ title = {{CDCS}: Cone Decomposition Conic Solver, version 1.1},
howpublished = {\url{https://github.com/giofantuzzi/CDCS}},
month = Sep,
year = 2016
}
```
A selection of BibTex styles that support arXiv preprints can be found [here](http://arxiv.org/hypertex/bibstyles/).
+
+
## Contact us
To contact us about CDCS, suggest improvements and report bugs, email either [Giovanni Fantuzzi](mailto:gf910@ic.ac.uk?Subject=CDCS) or [Yang Zheng](mailto:yang.zheng@eng.ox.ac.uk?Subject=CDCS).
diff --git a/README.txt b/README.txt
index e3cdf6e..c795563 100644
--- a/README.txt
+++ b/README.txt
@@ -4,11 +4,15 @@
CDCS (Cone Decomposition Conic Solver) is an open-source MATLAB solver for sparse
conic programs with partially decomposable conic constraints. CDCS implements the
-alternating direction method of multipliers (ADMM) described in on our paper:
-"Fast ADMM for Semidefinite Programs with Chordal Sparsity" available from
-(https://arxiv.org/pdf/1609.06068v2.pdf).
+alternating direction method of multipliers (ADMM) described in on our papers
+
+* "Fast ADMM for Semidefinite Programs with Chordal Sparsity" available from
+(https://arxiv.org/pdf/1609.06068v2.pdf)
-Current version: 1.0
+* "Fast ADMM for homogeneous self-dual embeddings of sparse SDPs" available from
+(https://arxiv.org/abs/1611.01828)
+
+Current version: 1.1
Release notes: CDCS is based on a temporary research code called ADMM-PDCP, which
is no longer maintained. If you downloaded ADMM-PDCP, please
replace it with CDCS.
@@ -32,7 +36,7 @@ Release notes: CDCS is based on a temporary research code called ADMM-PDCP, whic
CDCS solves in the standard primal and dual vectorized forms
minimize c'x maximize b'y
- (1) subject to Ax = b, (2) subject to c - A'y = z,
+ (1) subject to Ax = b, (2) subject to A'y + z = c,
x \in K z \in K*
where the conic constraint `x \in K` are partially decomposable. This means that
@@ -87,9 +91,17 @@ us (see the Contact Us section below).
HOW TO CITE
================================================================================
-If you find CDCS useful, please cite it with
+If you find CDCS useful, please cite
+
+@article{{ZFPGWhsde2016,
+ archivePrefix = {arXiv},
+ eprint = {1611.01828},
+ primaryClass = "math-OC",
+ author = {Zheng, Yang and Fantuzzi, Giovanni and Papachristodoulou, Antonis and Goulart, Paul and Wynn, Andrew},
+ title = {{Fast ADMM for homogeneous self-dual embeddings of sparse SDPs}}
+ }
-@article{ZFPGW2016,
+@article{ZFPGWpd2016,
archivePrefix = {arXiv},
eprint = {1609.06068v2},
primaryClass = "math-OC",
diff --git a/cdcs.m b/cdcs.m
index fe06407..eb3ca50 100644
--- a/cdcs.m
+++ b/cdcs.m
@@ -1,44 +1,50 @@
function [x,y,z,info] = cdcs(At,b,c,K,userOpts,initVars)
-
% CDCS
%
% Syntax:
%
-% [x,y,z,info] = CDCS(At,b,c,K,opts)
+% [x,y,z,info] = CDCS(At,b,c,K,options)
%
-% Solve a sparse conic program using chordal decomposition for the semidefinite
+% Solve a sparse conic program using chordal decomposition for the positive semidefinite
% cones and ADMM. CDCS solves the primal (P) or dual (D) standard forms of
% the conic problem,
%
% min max
-% (P) s.t. Ax = b, (D) s.t. c - A^Ty = z
+% (P) s.t. Ax = b, (D) s.t. A^Ty + z = c
% x \in K z \in K*
%
-% where A,b and c are the problem date and K is the cone (K* is the dual cone).
+% where A,b and c are the problem data and K is the cone (K* is the dual cone).
+% CDCS supports the following cones: Free, Linear, second-order,
+% Semi-definite, called as called K.f, K.l, K.q, and K.s.
+%
% The standard form to be solved is specified by the "solver" field of the
% options structure:
%
-% opts.solver = 'primal' (default): solve the problem in primal standard form
-% opts.solver = 'dual' : solve the problem in dual standard form
+% options.solver = 'hsde' (default): solve the problem in homogeneous self-dual embedding form
+% options.solver = 'primal' : solve the problem in primal standard form
+% options.solver = 'dual' : solve the problem in dual standard form
%
% The chordal decomposition can be carried out in two ways, specified by the
% "chordalize" option:
%
-% opts.chordalize = 1 (default): split the data equally between the cliques
-% opts.chordalize = 2 : assign data to one clique only
+% options.chordalize = 1 (default): split the data equally between the cliques
+% options.chordalize = 2 : assign data to one clique only
%
% Click here for a complete list of options.
%
% The output structure 'info' contains the following information:
%
% info.problem: - 0: CDCS terminated succesfully
-% - 1: the maximum number of iterations was reached
-% - 2: the ADMM iterations terminated succesfully, but the positive
-% matrix completion algorithm threw an error
+% - 1: primal infeasibility detected
+% - 2: dual infeasibility detected
+% - 3: maximum number of iterations reached
+% - 4: the ADMM iterations terminated succesfully, but the positive
+% matrix completion algorithm threw an error
% info.iter: number of iterations
% info.cost: terminal cost
% info.pres: terminal primal ADMM residual
% info.dres: terminal dual ADMM residual
+% info.log : history log of the ADMM residuals, cost, etc.
% info.time: some timing information (setup, ADMM iterations, cleanup, total)
%
% See also CDCSOPTS
@@ -61,10 +67,11 @@
%============================================
-% Solver options
+% Solver options & import cdcs_utils
%============================================
tstart = tic;
opts = cdcsOpts;
+import cdcs_utils.*
%============================================
@@ -76,13 +83,13 @@
end
% Checks on specified solver type and method
-if ~any(strcmpi(opts.solver,{'primal','dual'}))
- error('Unknown opts.solver. Please use "primal" or "dual".')
+if ~any(strcmpi(opts.solver,{'primal','dual','hsde'}))
+ error('Unknown opts.solver. Please use "primal", "dual" or "hsde".')
end
% Print nice welcoming header
if opts.verbose
- myline1 = [repmat('=',1,64),'\n'];
+ [header,myline1,myline2] = printHeader(opts);
fprintf(myline1)
fprintf('CDCS by G. Fantuzzi, Y. Zheng -- v1.0\n')
fprintf(myline1)
@@ -93,15 +100,14 @@
% start timing
proctime = tic;
-% sparsify everything, check cone constraints, rescale
+% sparsify everything, check cone constraints
[At,b,c,K,opts] = checkInputs(At,b,c,K,opts);
-[At,b,c,K,opts] = rescaleData(At,b,c,K,opts);
[At,b,c,K,opts] = splitBlocks(At,b,c,K,opts);
[opts.n,opts.m] = size(At);
-% chordal decomposition
+% rescale & chordal decomposition
Kold = K;
-[At,b,c,K,Ech,cd,chstuff] = chordalize(At,b,c,K,opts);
+[At,b,c,K,Ech,chstuff,opts] = preprocess(At,b,c,K,opts);
% basic decomposed problem dimensions: no. of cones, no. of vectorized conic
% variables, and no. of free primal variables
@@ -114,73 +120,69 @@
[X,Y,Z,others] = makeVariables(K,initVars,opts);
% Make operators for ADMM
-[step1,step2,step3,checkConv] = makeADMM(At,b,c,K,cd,Ech,opts);
+[updateX,updateY,updateZ,checkConvergence] = makeADMM(At,b,c,K,Ech,opts);
% Time setup and display
proctime = toc(proctime);
if opts.verbose
- myline2 = [repmat('-',1,64),'\n'];
+ % Set method to display
+ if strcmpi(opts.solver,'hsde')
+ method = 'homogeneous self-dual embedding';
+ else
+ method = opts.solver;
+ end
fprintf('done in %.4f seconds. \n',proctime);
- fprintf('Standard form : %s\n',opts.solver);
+ fprintf('Algorithm : %s\n',method);
fprintf('Chordalization method : %i\n',opts.chordalize);
fprintf('Adaptive penalty : %i\n',opts.adaptive);
fprintf('Scale data : %i\n',opts.rescale);
fprintf('Free variables : %i \n',K.f);
fprintf('Non-negative variables : %i \n',K.l);
- fprintf('Second-order cones : %i (max. size: %i)\n',length(K.q),max(K.q));
- fprintf('Semidefinite cones : %i (max. size: %i)\n',length(K.s),max(K.s));
+ fprintf('Second-order cones : %i (max. size: %i)\n',length(find(K.q ~=0)),max(K.q));
+ fprintf('Semidefinite cones : %i (max. size: %i)\n',length(find(K.s ~=0)),max(K.s));
fprintf('Affine constraints : %i \n',opts.m);
fprintf('Consensus constraints : %i \n',sum(accumarray(Ech,1)));
- fprintf(myline1)
+ fprintf(myline1);
+ fprintf(header);
+ fprintf(myline2);
end
%============================================
% Run ADMM
%============================================
-% Display
-if opts.verbose
- fprintf(' iter | pres | dres | cost | rho | time (s) |\n')
- fprintf(myline2)
-end
-
admmtime = tic;
-opts.feasCode = 1;
for iter = 1:opts.maxIter
% Save current iterate for convergence test
YOld = Y;
% Update block variables
- [X,others] = step1(X,Y,Z,opts.rho,others);
- [Y,others] = step2(X,Y,Z,opts.rho,others);
- [Z,others] = step3(X,Y,Z,opts.rho,others);
+ [X,others] = updateX(X,Y,Z,opts.rho,others);
+ [Y,others] = updateY(X,Y,Z,opts.rho,others);
+ [Z,others] = updateZ(X,Y,Z,opts.rho,others);
% log errors / check for convergence
- [isConverged,log,opts] = checkConv(X,Y,Z,YOld,others,iter,admmtime,opts);
- if isConverged
- opts.feasCode = 0;
+ [stop,info,log,opts] = checkConvergence(X,Y,Z,YOld,others,iter,admmtime,opts);
+ if stop
break;
end
end
admmtime = toc(admmtime);
-if opts.verbose
- fprintf(myline1)
-end
%============================================
% Outputs
%============================================
% Variables in sedumi format
posttime = tic;
-[x,y,z,opts] = setOutputs(X,Y,Z,others,Kold,cd,Ech,chstuff,opts);
+[x,y,z,info,opts] = setOutputs(X,Y,Z,others,Kold,c,Ech,chstuff,info,opts);
posttime = toc(posttime);
% Info
-info.problem = opts.feasCode; % diagnostic code
info.iter = iter; % # of iterations
-info.cost = log.cost; % terminal cost
-info.pres = log.pres; % terminal primal ADMM res
-info.dres = log.dres; % terminal dual ADMM res
+info.cost = log(iter).cost; % terminal cost
+info.pres = log(iter).pres; % terminal primal ADMM res
+info.dres = log(iter).dres; % terminal dual ADMM res
+info.log = log; % log of residuals etc
info.time.setup = proctime; % setup time
info.time.admm = admmtime; % ADMM time
info.time.cleanup = posttime; % post-processing time
@@ -188,6 +190,7 @@
% Print summary
if opts.verbose
+ fprintf(myline1)
fprintf(' SOLUTION SUMMARY:\n')
fprintf('------------------\n')
fprintf(' Termination code : %11.1d\n',info.problem)
diff --git a/cdcsInstall.m b/cdcsInstall.m
index bdfded1..7b30b0a 100644
--- a/cdcsInstall.m
+++ b/cdcsInstall.m
@@ -12,11 +12,13 @@
cd('include')
cs_install
cd(here);
-movefile(['include',filesep,'cs_lsolve.mex*'],[here,filesep,'private']);
-movefile(['include',filesep,'cs_ltsolve.mex*'],[here,filesep,'private']);
+movefile(['include',filesep,'cs_lsolve.mex*'], ...
+ [here,filesep,'packages',filesep,'+cdcs_utils']);
+movefile(['include',filesep,'cs_ltsolve.mex*'], ...
+ [here,filesep,'packages',filesep,'+cdcs_utils']);
% Then compile some mex files from this package
-cd('private')
+cd(['packages',filesep,'+cdcs_utils',filesep,'private'])
if (~isempty (strfind (computer, '64')))
mexcmd = 'mex -largeArrayDims' ;
else
@@ -27,6 +29,7 @@
cd(here)
% Finally add to path and save
+addpath([here,filesep,'packages']);
addpath(here);
savepath
diff --git a/cdcsOpts.m b/cdcsOpts.m
index 2eeb4ae..7298736 100644
--- a/cdcsOpts.m
+++ b/cdcsOpts.m
@@ -6,7 +6,7 @@
%
% Generic solver options
% ----------------------
-% options.solver = 'primal'; % which solver (primal/dual)?
+% options.solver = 'hsde'; % which solver (primal/dual/hsde)
% options.relTol = 1e-4; % tolerance
% options.rescale = true; % scale data to improve convergence
% options.verbose = 1; % print or silent
@@ -15,7 +15,7 @@
%
% Chordal decomposition options
% -----------------------------
-% options.chordalize = 1; % how to decompose the constraints (1/2)
+% options.chordalize = 1; % how to decompose the constraints (0/1/2)
% options.yPenalty = true; % add penalty term for Y block to cost
% options.completion = true; % complete the unused entries of the decomposed
% primal variable
@@ -43,7 +43,7 @@
% Create options structure
-options.solver = 'primal'; % which solver (primal/dual)?
+options.solver = 'hsde'; % which solver (primal/dual/hsde)
options.relTol = 1e-4; % tolerance
options.rescale = true; % scale data to improve convergence
options.verbose = 1; % print or silent
@@ -60,3 +60,7 @@
options.rhoMin = 1e-6; % minimum penalty parameter
options.rhoIt = 10; % if pres/dres>mu ( Y, v --> Z;
+% MAKEVARIABLES(K,initVars)
+% Make initial variables
+
+% GF:
+% Initialize with all zeros since we know it is a feasible point.
+% Initializing kappa=tau=1 gives annoying NaN in first iteration...
+
+
+if isempty(initVars)
+ %initialize vectorized versions of all variables
+ X.x = zeros(opts.nX,1)+1; % free variables
+ X.xh = zeros(opts.nXk,1)+1; % cone variables: vectorized version
+ X.y = zeros(opts.m,1); % free variables
+ X.v = zeros(opts.nXk,1); % free variables
+ X.tau = 0;
+
+ Y.x = zeros(opts.nX,1)+1; % free variables
+ Y.xh = zeros(opts.nXk,1)+1; % cone variables: vectorized version
+ Y.y = zeros(opts.m,1); % free variables
+ Y.v = zeros(opts.nXk,1); % free variables
+ Y.tau = 0;
+
+ Z.x = zeros(opts.nX,1); % similar to Y
+ Z.xh = zeros(opts.nXk,1); %
+ Z.y = zeros(opts.m,1); %
+ Z.v = zeros(opts.nXk,1); %
+ Z.kappa = 0;
+
+ others.X = cdcs_utils.makeConeVariables(K);% cone variables: matrix version
+
+else
+ % initialized vectorized variables with used input
+ % TO DO LATER
+
+end
+
+
+
+
+end
\ No newline at end of file
diff --git a/packages/+cdcs_hsde/preprocess.m b/packages/+cdcs_hsde/preprocess.m
new file mode 100644
index 0000000..80b8244
--- /dev/null
+++ b/packages/+cdcs_hsde/preprocess.m
@@ -0,0 +1,82 @@
+function [At,b,c,K,Ech,stuff,opts] = preprocess(At,b,c,K,opts)
+
+% CDCS/packages/+cdcs_hsde/PREPROCESS
+% Chordal decomposition of SDP constraints & rescale data for homogeneous
+% self-dual embedding
+
+% Import functions
+import cdcs_utils.svecData
+import cdcs_utils.chordalDecomposition
+
+% Any variables at all?
+nonSDP = K.f + K.l + sum(K.q);
+if nonSDP + sum(K.s) == 0
+ error('No variables in your problem?')
+end
+
+%--------------------------------------------
+% Chordal decomposition
+%--------------------------------------------
+
+% Returns the submatrix Ats, Cs of data for PSD variables with svec
+% instead of vec, which is the standard input. Also return modified
+% matrices At,c that account for svec operation.
+[At,c,Ats,Cs] = svecData(At,c,K);
+totvars = size(At,1);
+
+% If required and have SDP variables, decompose
+if ~isempty(K.s) && any(K.s~=0) && opts.chordalize~=0
+ [cliques,Ech,Jch,usedvars,s] = chordalDecomposition(Ats,Cs,K);
+
+else
+ cliques = [];
+ Ech = (1:totvars)';
+ Jch = (1:totvars)';
+ usedvars = (1:totvars)';
+ s = K.s;
+
+end
+
+%--------------------------------------------
+% Rescale data
+%--------------------------------------------
+[At,b,c,K,opts] = rescaleData(At,b,c,K,Jch,opts);
+
+%--------------------------------------------
+% Set new At, c, K.s (remove unused variables)
+%--------------------------------------------
+At = At(usedvars,:);
+c = c(usedvars);
+K.s = s;
+
+if opts.rescale == 1
+ opts.scaleFactors.D1 = opts.scaleFactors.D1(usedvars);
+end
+
+% Check if At,b,C are indeed sparse - if not, make full for speed!
+[n,m] = size(At);
+densityTol = 0.6;
+if nnz(At)/(m*n) > densityTol
+ At = full(At);
+end
+if nnz(b)/m > densityTol
+ b = full(b);
+end
+if nnz(c)/n > densityTol
+ c = full(c);
+end
+
+
+%--------------------------------------------
+% Set stuff
+%--------------------------------------------
+stuff.cliques = cliques;
+stuff.usedvars = usedvars;
+stuff.totvars = totvars;
+
+% norms, used repeatly in convergence checking
+opts.nAt = norm(At,'fro');
+opts.nb = norm(b,'fro');
+opts.nc = norm(c,'fro');
+
+
diff --git a/packages/+cdcs_hsde/printHeader.m b/packages/+cdcs_hsde/printHeader.m
new file mode 100644
index 0000000..30e7652
--- /dev/null
+++ b/packages/+cdcs_hsde/printHeader.m
@@ -0,0 +1,10 @@
+function [header,myline1,myline2] = printHeader
+
+% CDCS/packages/cdcs_pd/PRINTHEADER
+% Print header for iterations
+
+% Set stuff
+myline1 = [repmat('=',1,86),'\n'];
+myline2 = [repmat('-',1,86),'\n'];
+header = [' iter | pres | dres | pcost | dcost | gap ',...
+ '| rho | time (s) |\n'];
diff --git a/packages/+cdcs_hsde/private/checkConvergence.m b/packages/+cdcs_hsde/private/checkConvergence.m
new file mode 100644
index 0000000..6edd2b9
--- /dev/null
+++ b/packages/+cdcs_hsde/private/checkConvergence.m
@@ -0,0 +1,127 @@
+function [stop,info,log,opts] = checkConvergence(hatu,u,v,uold,iter,admmtime,opts,At,A,b,c,btr,ctr,E,others)
+
+% CHECKCONVERGENCE primal/dual infeasible
+% Problem codes
+%
+% 0: converged
+% 1: primal infeasible
+% 2: dual infeasible
+% 3: max number of iterations reached
+
+stop = false;
+info.problem = 3;
+
+if u.tau > 0 %% feasible problem
+ x = u.x./u.tau; y = u.y./u.tau;
+ dcost = (btr*y)/opts.scaleFactors.sc_cost; %% dual cost
+ pcost = (ctr*x)/opts.scaleFactors.sc_cost; %% primal cost
+
+ % this is based on the scaled data
+ presi = norm(A*x - b,'fro')./max([1,opts.nAt,opts.nb]);
+ dresi = norm((At*y + accumarray(E,u.v./u.tau)- c),'fro')./max([1,opts.nAt,opts.nc]);
+ gap = abs(pcost - dcost)/(1 + abs(pcost) + abs(dcost));
+
+ if max([presi,dresi,gap]) < opts.relTol
+ info.problem = 0;
+ stop = true;
+ end
+
+ % adpative penalty?
+ % GF: Only adapt penalty if not stopping after this iteration...
+ if opts.adaptive && ~stop
+
+ % standard ADMM residual according to Boyd's paper
+ utmp = vertcat(u.x,u.xh,u.y,u.v,u.tau);
+ vtmp = vertcat(v.x,v.xh,v.y,v.v,v.kappa);
+ hatutmp = vertcat(hatu.x,hatu.xh,hatu.y,hatu.v,hatu.tau);
+ uoldtmp = vertcat(uold.x,uold.xh,uold.y,uold.v,uold.tau);
+
+ r = norm(utmp-hatutmp,'fro')./max(norm(utmp,'fro'),norm(hatutmp,'fro'));
+ s = opts.rho*norm(utmp-uoldtmp,'fro')./max(norm(utmp,'fro'),norm(vtmp,'fro'));
+
+ % Update ADMM penalty parameter
+ opts = updatePenalty(opts,r,s);
+
+ end
+
+else % infeasible or unbounded problem?
+ pinfIndex = -btr*u.y; % index of certificating primal infesibility
+ dinfIndex = ctr*u.x; % index of certificating dual infesibility
+ if dinfIndex < 0 % dual infeasible
+ pfeas = norm(A*u.x,'fro');
+ pfeasTol = -dinfIndex/opts.nc*opts.relTol;
+ if pfeas <= pfeasTol %% a point x that certificates dual infeasiblity
+ info.problem = 2;
+ stop = true;
+ end
+ end
+ if pinfIndex < 0 % primal infeasible
+ dfeas = norm(At*u.y+accumarray(E,u.v),'fro');
+ dfeasTol = -pinfIndex/opts.nb*opts.relTol;
+ if dfeas <= dfeasTol %% a point y that certificates primal infeasiblity
+ info.problem = 1;
+ stop = true;
+ end
+ end
+
+ % value
+ presi = NaN;
+ dresi = NaN;
+ pcost = NaN;
+ dcost = NaN;
+ gap = NaN;
+end
+
+%progress message
+if opts.verbose && (iter == 1 || ~mod(iter,opts.dispIter) || stop)
+ fprintf('%5d | %7.2e | %7.2e | %9.2e | %9.2e | %8.2e | %8.2e | %8.2e |\n',...
+ iter,presi,dresi,pcost,dcost,gap, opts.rho,toc(admmtime))
+end
+
+% log information
+log(iter).pres = presi;
+log(iter).dres = dresi;
+log(iter).cost = pcost; %% use primal cost
+log(iter).dcost = dcost;
+
+% end main
+end
+
+% ============================================================================ %
+% Nested functions
+% ============================================================================ %
+
+function opts = updatePenalty(opts,pres,dres)
+% Update penaly
+
+ persistent itPinf itDinf
+ if isempty(itPinf) || isempty(itDinf)
+ % initialize persistent iteration counters when entering for the first
+ % time (persistent variables are empty and not zero when declared)
+ itPinf = 0; % number of iterations for which pinf/dinf <= eta
+ itDinf = 0; % number of iterations for which pinf/dinf > eta
+ end
+
+ if opts.adaptive
+ resRat = pres/dres;
+ if resRat >= opts.mu
+ itPinf = itPinf+1;
+ itDinf = 0;
+ if itPinf >= opts.rhoIt
+ % ratio of pinf and dinf remained large for long => rescale rho
+ itPinf = 0;
+ opts.rho = min(opts.rho*opts.tau, opts.rhoMax);
+ end
+ elseif 1/resRat >= opts.mu
+ itDinf = itDinf+1;
+ itPinf = 0;
+ if itDinf >= opts.rhoIt
+ % ratio of pinf and dinf remained small for long => rescale rho
+ itDinf = 0;
+ opts.rho = max(opts.rho/opts.tau, opts.rhoMin);
+ end
+ end
+ end
+
+end
+
diff --git a/packages/+cdcs_hsde/private/factorMatrix.m b/packages/+cdcs_hsde/private/factorMatrix.m
new file mode 100644
index 0000000..67629c2
--- /dev/null
+++ b/packages/+cdcs_hsde/private/factorMatrix.m
@@ -0,0 +1,150 @@
+function [xi,solInner] = factorMatrix(At,b,c,E,flag)
+% generate projector onto affine constraints
+% Here, we do not consider the modification of H in the scaling process
+
+D = accumarray(E,1);
+P = 1./(1+D./2);
+
+% First factor the matrix (I+APA') where P = inv(I+1/2D) diagonal
+factors = factorKKT(P,At,flag);
+
+% Compute a useful vector
+z.x = c;
+z.xh = sparse(length(E),1);
+z.y = -b;
+z.v = sparse(length(E),1);
+xi = solveInner(factors,E,z);
+const = 1+c'*xi.x - b'*xi.y;
+xi.x = xi.x/const;
+xi.xh = xi.xh/const;
+xi.y = xi.y/const;
+xi.v = xi.v/const;
+
+% Set function handle
+solInner = @(v)solveInner(factors,E,v);
+
+end
+
+%----------------------
+% solveInner
+%----------------------
+function u = solveInner(factors,E,v)
+ % Solve system of form
+ %
+ % [ I M1] [u1] = [v1]
+ % [-M1' I] [u2] [v2]
+ %
+ % where u1 = [u11,u12], u2 = [u21,u22] etc. Input v is formatted so that
+ % v.x = v11
+ % v.xh = v12
+ % v.y = v21
+ % v.v = v22
+ %
+ % Output is formatted the same way
+
+ persistent useBuiltin
+
+ if(isempty(useBuiltin))
+ %default to look for CSparse code
+ useBuiltin = ~exist(['cs_ltsolve.' mexext],'file');
+ end
+
+ % Import functions
+ import cdcs_utils.cs_lsolve
+ import cdcs_utils.cs_ltsolve
+
+ if strcmpi(factors.flag,'blk')
+
+ % Cholesky factors and original blocks
+ R = factors.R;
+ A = factors.A;
+ At = factors.At;
+ P = factors.P;
+ s = factors.s;
+ si = factors.si;
+
+ % First find v0 to solve system like M*u2 = v0
+ v01 = accumarray(E,v.v) + (At*v.y) + v.x;
+ v02 = v.xh - v.v;
+
+ % Solve system M*u2=v0 using factors
+ z = P.*(v01 + accumarray(E,v02)./2);
+ d = A*z;
+ ds= full(d(s,1)); % permute
+ if(useBuiltin)
+ p = R'\((R\ds)); % Native matlab version (slow)
+ else
+ p = cs_ltsolve(R,cs_lsolve(R,ds)); % Csparse version (avoids transpose)
+ end
+ p = p(si,1); % permute back
+
+ % finish solve
+ u.x = z - P.*(At*p);
+ u.xh = (v02 + u.x(E))./2;
+ u.y = v.y - A*u.x;
+ u.v = v.v + u.xh - u.x(E);
+ else
+ error('Unknown method')
+ end
+
+
+end
+
+
+%----------------------
+% factorKKT
+%----------------------
+function factors = factorKKT(P,At,flag)
+
+ if(nargin < 2)
+ At = [];
+ flag = 'blk';
+ end
+
+ if(nargin < 3)
+ flag = 'blk';
+ end
+
+ %[n,~] = size(At);
+
+ if strcmpi(flag,'ldl')
+ % TO DO
+ % I don't think ldl decomposition is a good choice for hsde.
+ error('Unknown method')
+
+ elseif strcmpi(flag,'inv')
+ % TO DO
+ error('Unknown method')
+ else
+ [n,m] = size(At);
+ B = spdiags(sqrt(P),0,n,n)*At; % slow with sparse type with few nonzeros?
+ M = speye(m) + sparse(B'*B); % = I + A*P*A'
+ [R,p,s] = chol(M,'lower','vector');
+
+ if p==0
+ %For maximum efficiency in projection, store both
+ %the permutation s and its inverse permutation
+ factors.flag = 'blk';
+ factors.R = R;
+ factors.A = At';
+ factors.At = At;
+ factors.P = P;
+ factors.s = s;
+ tmp = 1:length(s);
+ factors.si(s) = tmp; %inverse permutation
+ else
+ error('Cholesky decomposition fails ... ')
+ % not definite - use LDL
+ % M = [spdiags(1./P,0,n,n), sparse(At)
+ % sparse(m,n), sparse(m,m)];
+ % [U,D,s] = ldl(M,'upper','vector');
+ % factors.flag = 'ldl';
+ % factors.L = U';
+ % factors.D = D;
+ % factors.s = s;
+ % tmp = 1:length(s);
+ % factors.si(s) = tmp; %inverse permutation
+ end
+ end
+
+end
diff --git a/packages/+cdcs_hsde/private/rescaleData.m b/packages/+cdcs_hsde/private/rescaleData.m
new file mode 100644
index 0000000..8e43ea9
--- /dev/null
+++ b/packages/+cdcs_hsde/private/rescaleData.m
@@ -0,0 +1,124 @@
+function [At,b,c,K,opts] = rescaleData(At,b,c,K,Ech,opts)
+
+% +cdcs_hsde/RESCALEDATA.m
+% Try to rescale data to get nicer convergence properties.
+% taking the matching matrix H into account
+% Assumes no zero columns in At, otherwise will fail.
+
+% Parameters
+[n,m] = size(At);
+min_scale = 1e-3;
+max_scale = 1e+3;
+minScaleRowAt = min_scale .*sqrt(m);
+maxScaleRowAt = max_scale .*sqrt(m);
+minScaleColAt = min_scale .*sqrt(n);
+maxScaleColAt = max_scale .*sqrt(n);
+
+if opts.rescale
+
+ % -----------------------------------------------------
+ % Scale rows of [At H']
+ % [0 -I]
+ % -----------------------------------------------------
+
+ % norm of rows of [At H']
+ D1 = full(sqrt( sum(At.*At,2) + accumarray(Ech,1) ));
+ count = 0;
+
+ % No need to take mean for K.f and K.l
+ if K.f>0
+ count = count + K.f;
+ end
+
+ if K.l>0
+ count = count + K.l;
+ end
+
+ % mean over quaratic cones
+ if sum(K.q)>0
+ for i = 1:length(K.q)
+ nvars = K.q(i);
+ D1(count+1:count+nvars) = mean(D1(count+1:count+nvars));
+ count = count + nvars;
+ end
+ end
+
+ % mean over SDP cones (with svec!)
+ if sum(K.s)>0
+ for i = 1:length(K.s)
+ nvars = K.s(i)*(K.s(i)+1)/2;
+ D1(count+1:count+nvars) = mean(D1(count+1:count+nvars));
+ count = count + nvars;
+ end
+ end
+
+ D1(D1>maxScaleRowAt) = maxScaleRowAt; % set upper bound
+ D1(D1maxScaleColAt) = maxScaleColAt; % set upper bound
+ E1(E1 densityTol
+ At = full(At);
+end
+if nnz(b)/m > densityTol
+ b = full(b);
+end
+if nnz(c)/n > densityTol
+ c = full(c);
+end
+
+
+%--------------------------------------------
+% Set stuff
+%--------------------------------------------
+stuff.cliques = cliques;
+stuff.usedvars = usedvars;
+stuff.totvars = totvars;
\ No newline at end of file
diff --git a/packages/+cdcs_pd/printHeader.m b/packages/+cdcs_pd/printHeader.m
new file mode 100644
index 0000000..d8e5294
--- /dev/null
+++ b/packages/+cdcs_pd/printHeader.m
@@ -0,0 +1,9 @@
+function [header,myline1,myline2] = printHeader
+
+% CDCS/packages/cdcs_pd/PRINTHEADER
+% Print header for iterations
+
+% Set stuff
+myline1 = [repmat('=',1,64),'\n'];
+myline2 = [repmat('-',1,64),'\n'];
+header = ' iter | pres | dres | cost | rho | time (s) |\n';
diff --git a/private/checkConvergence.m b/packages/+cdcs_pd/private/checkConvergence.m
similarity index 79%
rename from private/checkConvergence.m
rename to packages/+cdcs_pd/private/checkConvergence.m
index 5358cdc..f06d67e 100644
--- a/private/checkConvergence.m
+++ b/packages/+cdcs_pd/private/checkConvergence.m
@@ -1,7 +1,11 @@
-function [isConverged,log,opts] = checkConvergence(X,Y,Z,YOld,others,b,c,E,iter,opts,admmtime)
+function [stop,info,log,opts] = checkConvergence(X,Y,Z,YOld,others,b,c,E,iter,opts,admmtime)
% CHECKCONVERGENCE
% Use the basic convergence test in the Boyd survey paper
+% Primal/dual solver cannot detect infeasibility so error codes can only be
+% 0: problem successfully solved
+% 3: max number of iterations reached
+
persistent itPinf itDinf
if iter == 1
@@ -11,6 +15,8 @@
itDinf = 0; % number of iterations for which pinf/dinf > eta
end
+% Import functions
+import cdcs_utils.flatten
% Extract some variables
Ex = flatten(Y.vec,X.blk);
@@ -47,13 +53,15 @@
%stopping criteria
if(max(pres,dres) 0)
K.q = K.q(K.q~=0);
nConeVars = nConeVars + sum(K.q);
else
K.q = 0;
end
-if isfield(K,'s')
+if (isfield(K,'s') && max(K.s) > 0)
K.s = K.s(K.s~=0);
nConeVars = nConeVars + sum(K.s.^2);
else
diff --git a/packages/+cdcs_utils/chordalDecomposition.m b/packages/+cdcs_utils/chordalDecomposition.m
new file mode 100644
index 0000000..00d9784
--- /dev/null
+++ b/packages/+cdcs_utils/chordalDecomposition.m
@@ -0,0 +1,92 @@
+function [cliques,Ech,Jch,usedvars,s] = chordalDecomposition(Ats,Cs,K)
+
+
+% CDCS/packages/+cdcs_utils/CHORDALDECOMPOSITION.m
+%
+% Compute chordal decomposition of SDP cones.
+
+%--------------------------------------------
+% Non SDP variables
+%--------------------------------------------
+nonSDP = K.f + K.l + sum(K.q);
+nonSDPind = (1:nonSDP)';
+
+
+%--------------------------------------------
+% Sparsity pattern, choral extension and maximal cliques
+%--------------------------------------------
+SP = spones(spones(Cs) + sparse(sum(spones(Ats),2))); % vector of 1s and 0s
+
+% Decompose each block
+nCones = length(K.s);
+count = 0;
+countE = 0;
+usedvars = [];
+s = cell(nCones,1);
+cliques = cell(nCones,1);
+Ech = {};
+Jch = {};
+for i = 1:nCones
+
+ % make block
+ % n = K.s(i)^2;
+ % Blk = reshape(SP(count+1:count+n),K.s(i),K.s(i));
+ n = K.s(i)*(K.s(i)+1)/2;
+ Blk = smat(SP(count+1:count+n));
+
+ % Chordal decomposition
+ cliques{i} = cliquesFromSpMatD(Blk);
+ P = cliques{i}.NoC;
+ Blk = spones(cliques{i}.idxMatrix+cliques{i}.idxMatrix'); %% chordal sparsity
+ ri = cliques{i}.Elem;
+ ci = zeros(length(ri),1);
+ tmp = 0;
+ for k = 1:P
+ ci( tmp+1 : tmp+cliques{i}.NoElem(k) ) = k;
+ tmp = tmp + cliques{i}.NoElem(k);
+ end
+ MCO = sparse(ri,ci,1,K.s(i),P);
+
+ % Find variables to keep
+ nnzind = find(svec(Blk));
+ usedvars = [usedvars; nnzind+nonSDP+count];
+
+
+ % Indexing to extract local submatrices & split cone
+ E = cell(P,1);
+ ind = cell(P,1);
+ s{i} = zeros(1,P);
+ for j = 1:P
+
+ %indices in global variable X (including zero entries)
+ Position = find(MCO(:,j));
+ s{i}(j) = length(Position); % size of clique cone
+ Position = repmat(Position,1,s{i}(j));
+ rw = Position(:);
+ cl = Position'; cl = cl(:);
+ lti = rw>=cl; % lower triangular indices
+ rw = rw(lti);
+ cl = cl(lti);
+ ind{j} = rw+(K.s(i)-cl./2).*(cl-1); % linear indices in svec(X)
+
+ % Which entry of list of nonzero elements in global variable?
+ [LIA,LOCB] = ismember(ind{j},nnzind);
+ E{j} = LOCB + nonSDP + countE;
+ ind{j} = ind{j} + nonSDP + count;
+ end
+
+ Ech = [Ech; E];
+ Jch = [Jch; ind];
+ count = count + n;
+ countE = countE + length(nnzind);
+end
+
+
+
+%--------------------------------------------
+% Concatenate mapping indices
+%--------------------------------------------
+Ech = [nonSDPind; vertcat(Ech{:})];
+Jch = [nonSDPind; vertcat(Jch{:})];
+s = horzcat(s{:});
+usedvars = [nonSDPind; usedvars];
\ No newline at end of file
diff --git a/private/clean.m b/packages/+cdcs_utils/clean.m
similarity index 100%
rename from private/clean.m
rename to packages/+cdcs_utils/clean.m
diff --git a/private/cs_lsolve.m b/packages/+cdcs_utils/cs_lsolve.m
similarity index 100%
rename from private/cs_lsolve.m
rename to packages/+cdcs_utils/cs_lsolve.m
diff --git a/private/cs_ltsolve.m b/packages/+cdcs_utils/cs_ltsolve.m
similarity index 100%
rename from private/cs_ltsolve.m
rename to packages/+cdcs_utils/cs_ltsolve.m
diff --git a/private/flatten.m b/packages/+cdcs_utils/flatten.m
similarity index 100%
rename from private/flatten.m
rename to packages/+cdcs_utils/flatten.m
diff --git a/packages/+cdcs_utils/makeADMM.m b/packages/+cdcs_utils/makeADMM.m
new file mode 100644
index 0000000..2277bfc
--- /dev/null
+++ b/packages/+cdcs_utils/makeADMM.m
@@ -0,0 +1,22 @@
+function [updateX,updateY,updateZ,checkConvergence] = makeADMM(At,b,c,K,Ech,opts)
+
+% CDCS/packages/+cdcs_utils/MAKEADMM.m
+%
+% Initialize ADMM operators for CDCS. Call the correct initialization routine
+% according to the method specified in opts.solver.
+
+switch lower(opts.solver)
+
+ case {'primal', 'dual'}
+ [updateX,updateY,updateZ,checkConvergence] = ...
+ cdcs_pd.makeADMM(At,b,c,K,Ech,opts);
+
+ case {'hsde'}
+ [updateX,updateY,updateZ,checkConvergence] = ...
+ cdcs_hsde.makeADMM(At,b,c,K,Ech,opts);
+
+ otherwise
+ error('Unknown value for ''options.solver''.')
+
+end
+
diff --git a/private/makeConeVariables.m b/packages/+cdcs_utils/makeConeVariables.m
similarity index 100%
rename from private/makeConeVariables.m
rename to packages/+cdcs_utils/makeConeVariables.m
diff --git a/packages/+cdcs_utils/makeVariables.m b/packages/+cdcs_utils/makeVariables.m
new file mode 100644
index 0000000..3072504
--- /dev/null
+++ b/packages/+cdcs_utils/makeVariables.m
@@ -0,0 +1,20 @@
+function [X,Y,Z,others] = makeVariables(K,initVars,opts)
+
+% CDCS/packages/+cdcs_utils/MAKEVARIABLES.m
+%
+% Initialize variables for CDCS. Call the correct variable initialization
+% routine according to the method specified in opts.solver.
+
+switch lower(opts.solver)
+
+ case {'primal', 'dual'}
+ [X,Y,Z,others] = cdcs_pd.makeVariables(K,initVars,opts);
+
+ case {'hsde'}
+ %error('Homogeneous self-dual embedding solver coming soon!')
+ [X,Y,Z,others] = cdcs_hsde.makeVariables(K,initVars,opts);
+
+ otherwise
+ error('Unknown value for ''options.solver''.')
+
+end
\ No newline at end of file
diff --git a/packages/+cdcs_utils/preprocess.m b/packages/+cdcs_utils/preprocess.m
new file mode 100644
index 0000000..32b88c8
--- /dev/null
+++ b/packages/+cdcs_utils/preprocess.m
@@ -0,0 +1,20 @@
+function [At,b,c,K,Ech,chstuff,opts] = preprocess(At,b,c,K,opts)
+
+% CDCS/packages/+cdcs_utils/PREPROCESS.m
+%
+% Preprocess data: chordalize and rescale. Call functions according to their
+% method (homogeneous self-dual method uses different rescaling!)
+
+switch lower(opts.solver)
+
+ case {'primal', 'dual'}
+ [At,b,c,K,Ech,chstuff,opts] = cdcs_pd.preprocess(At,b,c,K,opts);
+
+ case {'hsde'}
+ [At,b,c,K,Ech,chstuff,opts] = cdcs_hsde.preprocess(At,b,c,K,opts);
+
+ otherwise
+ error('Unknown value for ''options.solver''.')
+
+end
+
diff --git a/packages/+cdcs_utils/printHeader.m b/packages/+cdcs_utils/printHeader.m
new file mode 100644
index 0000000..d09f07b
--- /dev/null
+++ b/packages/+cdcs_utils/printHeader.m
@@ -0,0 +1,20 @@
+function [header,myline1,myline2] = printHeader(opts)
+
+% CDCS/packages/+cdcs_utils/PRINTHEADER.m
+%
+% Print header for solver
+
+switch lower(opts.solver)
+
+ case {'primal', 'dual'}
+ [header,myline1,myline2] = cdcs_pd.printHeader;
+
+ case {'hsde'}
+ %error('Homogeneous self-dual embedding solver coming soon!')
+ [header,myline1,myline2] = cdcs_hsde.printHeader;
+
+ otherwise
+ error('Unknown value for ''options.solver''.')
+
+end
+
diff --git a/private/cliquesFromSpMatD.m b/packages/+cdcs_utils/private/cliquesFromSpMatD.m
similarity index 100%
rename from private/cliquesFromSpMatD.m
rename to packages/+cdcs_utils/private/cliquesFromSpMatD.m
diff --git a/private/constructCliqueGraph.m b/packages/+cdcs_utils/private/constructCliqueGraph.m
similarity index 100%
rename from private/constructCliqueGraph.m
rename to packages/+cdcs_utils/private/constructCliqueGraph.m
diff --git a/private/maxSpanningTree.m b/packages/+cdcs_utils/private/maxSpanningTree.m
similarity index 100%
rename from private/maxSpanningTree.m
rename to packages/+cdcs_utils/private/maxSpanningTree.m
diff --git a/private/smat.c b/packages/+cdcs_utils/private/smat.c
similarity index 100%
rename from private/smat.c
rename to packages/+cdcs_utils/private/smat.c
diff --git a/private/svec.c b/packages/+cdcs_utils/private/svec.c
similarity index 100%
rename from private/svec.c
rename to packages/+cdcs_utils/private/svec.c
diff --git a/private/projectK.m b/packages/+cdcs_utils/projectK.m
similarity index 100%
rename from private/projectK.m
rename to packages/+cdcs_utils/projectK.m
diff --git a/private/psdCompletion.m b/packages/+cdcs_utils/psdCompletion.m
similarity index 100%
rename from private/psdCompletion.m
rename to packages/+cdcs_utils/psdCompletion.m
diff --git a/packages/+cdcs_utils/setOutputs.m b/packages/+cdcs_utils/setOutputs.m
new file mode 100644
index 0000000..0fce28b
--- /dev/null
+++ b/packages/+cdcs_utils/setOutputs.m
@@ -0,0 +1,20 @@
+function [x,y,z,info,opts] = setOutputs(X,Y,Z,others,Kold,c,Ech,chstuff,info,opts)
+
+% CDCS/packages/+cdcs_utils/SETOUTPUTS.m
+%
+% Set output for CDCS. Call the correct routine according to the method
+% specified in opts.solver.
+
+switch lower(opts.solver)
+
+ case {'primal', 'dual'}
+ [x,y,z,info,opts] = cdcs_pd.setOutputs(X,Y,Z,others,Kold,c,Ech,chstuff,info,opts);
+
+ case {'hsde'}
+ %error('Homogeneous self-dual embedding solver coming soon!')
+ [x,y,z,info,opts] = cdcs_hsde.setOutputs(X,Y,Z,others,Kold,c,Ech,chstuff,info,opts);
+
+ otherwise
+ error('Unknown value for ''options.solver''.')
+
+end
\ No newline at end of file
diff --git a/private/setUserOpts.m b/packages/+cdcs_utils/setUserOpts.m
similarity index 87%
rename from private/setUserOpts.m
rename to packages/+cdcs_utils/setUserOpts.m
index 0547dfd..d1363a6 100644
--- a/private/setUserOpts.m
+++ b/packages/+cdcs_utils/setUserOpts.m
@@ -4,7 +4,6 @@
%
% Set user options
-
if ~isstruct(userOpts)
error('Input ''options'' must be a structure.');
else
@@ -16,4 +15,7 @@
warning('Option ''%s'' is unknown and will be ignored.',fnames{n})
end
end
-end
\ No newline at end of file
+end
+
+% make solver lowercase
+opts.solver = lower(opts.solver);
\ No newline at end of file
diff --git a/private/splitBlocks.m b/packages/+cdcs_utils/splitBlocks.m
similarity index 100%
rename from private/splitBlocks.m
rename to packages/+cdcs_utils/splitBlocks.m
diff --git a/packages/+cdcs_utils/svecData.m b/packages/+cdcs_utils/svecData.m
new file mode 100644
index 0000000..31a9942
--- /dev/null
+++ b/packages/+cdcs_utils/svecData.m
@@ -0,0 +1,75 @@
+function [At,C,Ats,Cs] = svecData(At,C,K)
+
+% SVECDATA.m
+%
+% Transform data from vectorize format to svectorized format for semidefinite
+% variables according to the cone K.
+
+%--------------------------------------------
+% Some variables
+%--------------------------------------------
+nSDP = sum(K.s.^2);
+nonSDP = K.f + K.l + sum(K.q);
+
+if nSDP > 0
+ %--------------------------------------------
+ % Build matrix Q such that svec(X)=Q*vec(X).
+ % Consider each block separately.
+ %--------------------------------------------
+ nBlk = length(K.s);
+ isr2 = 1./sqrt(2);
+ I = zeros(nSDP,1);
+ J = zeros(nSDP,1);
+ V = zeros(nSDP,1);
+ count = 0;
+ bdrow = 0;
+ bdcol = 0;
+
+ % Loop over blocks
+ for blk = 1:nBlk
+ n = K.s(blk);
+ rowcnt = 0;
+
+ % Loop over columns
+ for j=1:n
+
+ % Find indices
+ nels = n-j+1;
+ rws = repmat(1:nels,2,1);
+ rws = rws(:);
+ cls = [n*(j-1)+(j+1:n); n*((j+1:n)-1)+j];
+ I(count+1:count+2*nels-1) = rowcnt + rws(2:end) + bdrow;
+ J(count+1:count+2*nels-1) = [n*(j-1)+j; cls(:)] + bdcol;
+ V(count+1:count+2*nels-1) = [1; repmat(isr2,2*nels-2,1)];
+
+ % update counters
+ rowcnt = rowcnt + nels;
+ count = count + 2*nels-1;
+ end
+
+ % Update block diagonal counter
+ bdrow = bdrow + n*(n+1)/2;
+ bdcol = bdcol + n^2;
+ end
+
+ % Build sparse Q
+ Q = sparse(I,J,V,bdrow,nSDP);
+
+ % svec the data
+ Ats = Q*At(nonSDP+1:end,:);
+ Cs = Q*C(nonSDP+1:end,:);
+
+else
+ Ats = [];
+ Cs = [];
+
+end
+
+%--------------------------------------------
+% Modify At and c
+%--------------------------------------------
+At = [At(1:nonSDP,:); Ats];
+C = [ C(1:nonSDP); Cs];
+
+
+end
\ No newline at end of file
diff --git a/private/vec.m b/packages/+cdcs_utils/vec.m
similarity index 100%
rename from private/vec.m
rename to packages/+cdcs_utils/vec.m
diff --git a/private/chordalize.m b/private/chordalize.m
deleted file mode 100644
index bb1c893..0000000
--- a/private/chordalize.m
+++ /dev/null
@@ -1,159 +0,0 @@
-function [At,b,c,K,Ech,cd,stuff] = chordalize(At,b,c,K,opts)
-
-% Chordal decomposition of SDP constraints
-
-% Any variables at all?
-if K.f + K.l + sum(K.q) + sum(K.s) == 0
- error('No variables in your problem?')
-end
-
-% Sort out variables not in SDP cones
-nonSDP = K.f + K.l + sum(K.q);
-nonSDPind = (1:nonSDP)';
-
-
-% If have SDP variables to decompose
-if ~isempty(K.s) && any(K.s~=0)
-
- %--------------------------------------------
- % Separate part of At, C for PSD vars
- %--------------------------------------------
- % Returns the submatrix Ats, Cs of data for PSD variables with svec
- % instead of vec, which is the standard input. Also return modified
- % matrices At,c that account for svec operation.
- [At,c,Ats,Cs] = svecData(At,c,K);
- totvars = size(At,1);
-
- %--------------------------------------------
- % Sparsity pattern, choral extension and maximal cliques
- %--------------------------------------------
- SP = spones(spones(Cs) + sparse(sum(spones(Ats),2))); % vector of 1s and 0s
-
- % Decompose each block
- nCones = length(K.s);
- count = 0;
- countE = 0;
- usedvars = [];
- s = cell(nCones,1);
- cliques = cell(nCones,1);
- Ech = {};
- for i = 1:nCones
-
- % make block
-% n = K.s(i)^2;
-% Blk = reshape(SP(count+1:count+n),K.s(i),K.s(i));
- n = K.s(i)*(K.s(i)+1)/2;
- Blk = smat(SP(count+1:count+n));
-
- % Chordal decomposition
- cliques{i} = cliquesFromSpMatD(Blk);
- P = cliques{i}.NoC;
- Blk = spones(cliques{i}.idxMatrix+cliques{i}.idxMatrix'); %% chordal sparsity
- ri = cliques{i}.Elem;
- ci = zeros(length(ri),1);
- tmp = 0;
- for k = 1:P
- ci( tmp+1 : tmp+cliques{i}.NoElem(k) ) = k;
- tmp = tmp + cliques{i}.NoElem(k);
- end
- MCO = sparse(ri,ci,1,K.s(i),P);
-
- % Find variables to keep
- nnzind = find(svec(Blk));
- usedvars = [usedvars; nnzind+nonSDP+count];
-
-
- % Indexing to extract local submatrices & split cone
- E = cell(P,1);
- ind = cell(P,1);
- s{i} = zeros(1,P);
- for j = 1:P
-
- %indices in global variable X (including zero entries)
- Position = find(MCO(:,j));
- s{i}(j) = length(Position); % size of clique cone
- Position = repmat(Position,1,s{i}(j));
- rw = Position(:);
- cl = Position'; cl = cl(:);
- lti = rw>=cl; % lower triangular indices
- rw = rw(lti);
- cl = cl(lti);
- ind{j} = rw+(K.s(i)-cl./2).*(cl-1); % linear indices in svec(X)
-
- % Which entry of list of nonzero elements in global variable?
- [LIA,LOCB] = ismember(ind{j},nnzind);
- E{j} = LOCB + nonSDP + countE;
-
- end
-
- Ech = [Ech; E];
- count = count + n;
- countE = countE + length(nnzind);
- end
-
-
- %--------------------------------------------
- % Remove rows from At, C corresponding to variables that can be dropped
- % according to the aggregate sparsity pattern in SP
- %--------------------------------------------
- K.s = horzcat(s{:});
- usedvars = [(1:nonSDP)';usedvars];
- At = At(usedvars,:);
- c = c(usedvars);
-
- % Check if At,b,C are indeed sparse - if not, make full for speed!
- [n,m]=size(At);
- densityTol = 0.6;
- if nnz(At)/(m*n) > densityTol
- At = full(At);
- end
- if nnz(b)/m > densityTol
- b = full(b);
- end
- if nnz(c)/n > densityTol
- c = full(c);
- end
-
- %--------------------------------------------
- % Concatenate mapping indices
- %--------------------------------------------
- Ech = [nonSDPind; vertcat(Ech{:})];
-
- %--------------------------------------------
- % Decompose problem data if necessary
- %--------------------------------------------
-
- if opts.chordalize == 1
- % Decompose equally
- IA = accumarray(Ech,1);
- cd = c./IA; cd = cd(Ech);
-
- elseif opts.chordalize == 2
- % Decompose using only last entry
- nv = length(Ech);
- cd = zeros(nv,1);
- [U,IA] = unique(Ech,'last');
- cd(IA,:) = c(U,:);
-
- else
- error('Unknown chordal decomposition method.')
-
- end
-
-
-else
- Ech = nonSDPind;
- cliques = [];
- cd = c;
- usedvars = nonSDPind;
- totvars = size(At,1);
-
-end
-
-%--------------------------------------------
-% Set stuff
-%--------------------------------------------
-stuff.cliques = cliques;
-stuff.usedvars = usedvars;
-stuff.totvars = totvars;
-end
\ No newline at end of file
diff --git a/private/makeADMM.m b/private/makeADMM.m
deleted file mode 100644
index 55140df..0000000
--- a/private/makeADMM.m
+++ /dev/null
@@ -1,22 +0,0 @@
-function [step1,step2,step3,checkConv] = makeADMM(At,b,c,K,cd,Ech,opts)
-
-% Make ADMM operators
-
-
-% Useful stuff
-[projAffine,projCone] = makeProjectors(At,b,c,K,cd,Ech,opts);
-
-% Update steps
-step1 = @(X,Y,Z,rho,others)updateX(X,Y,Z,rho,others,Ech,K,projAffine,opts);
-step2 = @(X,Y,Z,rho,others)updateY(X,Y,Z,rho,others,projCone);
-step3 = @(X,Y,Z,rho,others)updateZ(X,Y,Z,rho,others,K);
-
-% Convergence check
-checkConv = @(X,Y,Z,YOld,others,iter,admmtime,opts)...
- checkConvergence(X,Y,Z,YOld,others,b,c,Ech,iter,opts,admmtime);
-
-
-
-end
-
-
diff --git a/private/svecData.m b/private/svecData.m
deleted file mode 100644
index 3ae6850..0000000
--- a/private/svecData.m
+++ /dev/null
@@ -1,67 +0,0 @@
-function [At,C,Ats,Cs] = svecData(At,C,K)
-
-% SVECDATA.m
-%
-% Transform data from vectorize format to svectorized format for semidefinite
-% variables according to the cone K.
-
-%--------------------------------------------
-% Some variables
-%--------------------------------------------
-n = size(At,1);
-nSDP = sum(K.s.^2);
-nonSDP = K.f + K.l + sum(K.q);
-
-%--------------------------------------------
-% Build matrix Q such that svec(X)=Q*vec(X).
-% Consider each block separately.
-%--------------------------------------------
-nBlk = length(K.s);
-isr2 = 1./sqrt(2);
-I = zeros(nSDP,1);
-J = zeros(nSDP,1);
-V = zeros(nSDP,1);
-count = 0;
-bdrow = 0;
-bdcol = 0;
-
-% Loop over blocks
-for blk = 1:nBlk
- n = K.s(blk);
- rowcnt = 0;
-
- % Loop over columns
- for j=1:n
-
- % Find indices
- nels = n-j+1;
- rws = repmat(1:nels,2,1);
- rws = rws(:);
- cls = [n*(j-1)+(j+1:n); n*((j+1:n)-1)+j];
- I(count+1:count+2*nels-1) = rowcnt + rws(2:end) + bdrow;
- J(count+1:count+2*nels-1) = [n*(j-1)+j; cls(:)] + bdcol;
- V(count+1:count+2*nels-1) = [1; repmat(isr2,2*nels-2,1)];
-
- % update counters
- rowcnt = rowcnt + nels;
- count = count + 2*nels-1;
- end
-
- % Update block diagonal counter
- bdrow = bdrow + n*(n+1)/2;
- bdcol = bdcol + n^2;
-end
-
-% Build sparse Q
-Q = sparse(I,J,V,bdrow,nSDP);
-
-%--------------------------------------------
-% Modify At and c
-%--------------------------------------------
-Ats = Q*At(nonSDP+1:end,:);
-Cs = Q*C(nonSDP+1:end,:);
-At = [At(1:nonSDP,:); Ats];
-C = [ C(1:nonSDP); Cs];
-
-
-end
\ No newline at end of file