[37] | 1 | function output = callmaxdet(interfacedata); |
---|
| 2 | |
---|
| 3 | % Author Johan Löfberg |
---|
| 4 | % $Id: callmaxdet.m,v 1.3 2006/07/25 19:02:39 joloef Exp $ |
---|
| 5 | |
---|
| 6 | [F_struc,F_blksz,G_struc,G_blksz] = sedumi2maxdet(interfacedata.F_struc,interfacedata.K); |
---|
| 7 | |
---|
| 8 | c = interfacedata.c; |
---|
| 9 | options = interfacedata.options; |
---|
| 10 | x0 = interfacedata.x0; |
---|
| 11 | onlyfeasible = (nnz(c)==0) & isempty(G_blksz); |
---|
| 12 | |
---|
| 13 | if isempty(G_blksz) |
---|
| 14 | G_blksz = 1; |
---|
| 15 | G_struc = [1 zeros(1,length(c))]; |
---|
| 16 | end |
---|
| 17 | |
---|
| 18 | if isempty(F_blksz) |
---|
| 19 | x = zeros(length(c),1)+NaN; |
---|
| 20 | problem = 11; |
---|
| 21 | infostr = yalmiperror(problems,'MAXDET failed since there are no inequality constraints in F(x)'); |
---|
| 22 | solveroutput = []; |
---|
| 23 | solvertime = 0; |
---|
| 24 | D_struc = []; |
---|
| 25 | end |
---|
| 26 | |
---|
| 27 | problem = 0; |
---|
| 28 | solvertimephase1=0; |
---|
| 29 | D_struc = []; |
---|
| 30 | if isempty(x0) |
---|
| 31 | solvertimephase1 = clock; |
---|
| 32 | showprogress('Calling MAXDET/phase1',options.showprogress); |
---|
| 33 | [x0,z0,w0,problem,infostr,solveroutput] = callmaxdetphase1(full(F_struc),full(F_blksz), full(G_struc),full(G_blksz), full(c), options); |
---|
| 34 | solvertimephase1 = etime(clock,solvertimephase1); |
---|
| 35 | end |
---|
| 36 | if (problem~=0) | (onlyfeasible==1) |
---|
| 37 | if isempty(x0) |
---|
| 38 | x = repmat(nan,length(c),1); |
---|
| 39 | else |
---|
| 40 | x = x0; |
---|
| 41 | end |
---|
| 42 | solvertime = solvertimephase1; |
---|
| 43 | else |
---|
| 44 | z0=zeros(size(F_struc,1),1); |
---|
| 45 | w0=zeros(size(G_struc,1),1); |
---|
| 46 | solvertime = clock; |
---|
| 47 | showprogress('Calling MAXDET',options.showprogress); |
---|
| 48 | if options.verbose==0 |
---|
| 49 | evalc('[x,Z,W,ul,hist,infostr]=maxdet(full(F_struc),F_blksz,full(G_struc), G_blksz,c'',x0,z0,w0,options.maxdet.AbsTol,options.maxdet.RelTol,options.maxdet.gam,options.maxdet.NTiters);'); |
---|
| 50 | else |
---|
| 51 | [x,Z,W,ul,hist,infostr]=maxdet(full(F_struc),F_blksz,full(G_struc), G_blksz,c',x0,z0,w0,options.maxdet.AbsTol,options.maxdet.RelTol,options.maxdet.gam,options.maxdet.NTiters); |
---|
| 52 | end |
---|
| 53 | solvertime = etime(clock,solvertime)+solvertimephase1; |
---|
| 54 | |
---|
| 55 | D_struc = [Z;W]; |
---|
| 56 | |
---|
| 57 | if options.savesolveroutput |
---|
| 58 | solveroutput.x = x; |
---|
| 59 | solveroutput.Z = Z; |
---|
| 60 | solveroutput.W = W; |
---|
| 61 | solveroutput.ul = ul; |
---|
| 62 | solveroutput.hist = hist; |
---|
| 63 | solveroutput.infostr = infostr; |
---|
| 64 | else |
---|
| 65 | solveroutput = []; |
---|
| 66 | end |
---|
| 67 | |
---|
| 68 | switch infostr |
---|
| 69 | case {'target reached','relative accuracy reached','absolute accuracy reached','absolute tolerance reached','relative tolerance reached'} |
---|
| 70 | problem = 0; |
---|
| 71 | case 'maximum Newton iteration exceeded' |
---|
| 72 | problem = 3; |
---|
| 73 | otherwise |
---|
| 74 | problem = 1; |
---|
| 75 | end |
---|
| 76 | infostr = yalmiperror(problem,'MAXDET'); |
---|
| 77 | end |
---|
| 78 | |
---|
| 79 | |
---|
| 80 | % Standard interface |
---|
| 81 | output.Primal = x; |
---|
| 82 | output.Dual = D_struc; |
---|
| 83 | output.Slack = []; |
---|
| 84 | output.problem = problem; |
---|
| 85 | output.infostr = infostr; |
---|
| 86 | output.solverinput = []; |
---|
| 87 | output.solveroutput= solveroutput; |
---|
| 88 | output.solvertime = solvertime; |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | |
---|
| 112 | |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | |
---|
| 116 | |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | |
---|
| 123 | function [x0,z0,w0,problem,infostr,solveroutput] = callmaxdetphase1(F_struc,F_blksz, G_struc,G_blksz, c, options); |
---|
| 124 | |
---|
| 125 | % Author Johan Löfberg |
---|
| 126 | % $Id: callmaxdet.m,v 1.3 2006/07/25 19:02:39 joloef Exp $ |
---|
| 127 | |
---|
| 128 | try |
---|
| 129 | if options.verbose==0 |
---|
| 130 | evalc('[x0,z0,w0,ul] = phase1(full(F_struc),F_blksz,full(G_struc), G_blksz,options.maxdet.gam,options.maxdet.AbsTol,options.maxdet.RelTol,options.maxdet.NTiters);'); |
---|
| 131 | else |
---|
| 132 | [x0,z0,w0,ul] = phase1(full(F_struc),F_blksz,full(G_struc), G_blksz,options.maxdet.gam,options.maxdet.AbsTol,options.maxdet.RelTol,options.maxdet.NTiters); |
---|
| 133 | end |
---|
| 134 | if ul(1)<0 |
---|
| 135 | problem = 0; |
---|
| 136 | elseif ul(2)>=0 |
---|
| 137 | problem = 1; |
---|
| 138 | else |
---|
| 139 | problem = 8; |
---|
| 140 | end |
---|
| 141 | % % Problems currently detected outside |
---|
| 142 | % problem = (ul(1)>0); |
---|
| 143 | |
---|
| 144 | % In case of problems, phase1 outputs the extended vector instead |
---|
| 145 | % of last primal iterate |
---|
| 146 | if problem |
---|
| 147 | x0 = x0(1:length(c)); |
---|
| 148 | end |
---|
| 149 | |
---|
| 150 | catch |
---|
| 151 | problem = 9; |
---|
| 152 | x0 = []; |
---|
| 153 | z0 = []; |
---|
| 154 | w0 = []; |
---|
| 155 | ul = []; |
---|
| 156 | end |
---|
| 157 | infostr = yalmiperror(problem,'MAXDET/phase1'); |
---|
| 158 | if options.savesolveroutput |
---|
| 159 | solveroutput.x = x0; |
---|
| 160 | solveroutput.Z = z0; |
---|
| 161 | solveroutput.W = w0; |
---|
| 162 | solveroutput.ul = ul; |
---|
| 163 | solveroutput.infostr = infostr; |
---|
| 164 | else |
---|
| 165 | solveroutput = []; |
---|
| 166 | end |
---|
| 167 | |
---|
| 168 | % CODE TAKEN DIRECTLY FROM MAXDET/PHASE1 |
---|
| 169 | % THE REASON FOR HAVING A COPY HERE IS THAT |
---|
| 170 | % A FILE NAMED PHASE1 EXIST FOR SP ALSO |
---|
| 171 | %MAXDET, version ALPHA, April 1996. |
---|
| 172 | |
---|
| 173 | %COPYRIGHT 1996 SHAO-PO WU, LIEVEN VANDENBERGHE AND STEPHEN BOYD |
---|
| 174 | %Permission to use, copy, modify, and distribute this software for |
---|
| 175 | %any purpose without fee is hereby granted, provided that this entire |
---|
| 176 | %notice is included in all copies of any software which is or includes |
---|
| 177 | %a copy or modification of this software and in all copies of the |
---|
| 178 | %supporting documentation for such software. |
---|
| 179 | %This software is being provided "as is", without any express or |
---|
| 180 | %implied warranty. In particular, the authors do not make any |
---|
| 181 | %representation or warranty of any kind concerning the merchantability |
---|
| 182 | %of this software or its fitness for any particular purpose. |
---|
| 183 | function [x,Z,W,ul] = phase1(F,F_blkszs,G,G_blkszs,gam,abstol,reltol,NTiters); |
---|
| 184 | |
---|
| 185 | m = size(F,2)-1; |
---|
| 186 | if (m ~= size(G,2)-1) |
---|
| 187 | error('F and G must have the same number of columns.'); |
---|
| 188 | end |
---|
| 189 | if (size(F,1) ~= sum(F_blkszs.^2)) |
---|
| 190 | error('Dimensions of F do not match F_blkszs.'); |
---|
| 191 | end; |
---|
| 192 | if (size(G,1) ~= sum(G_blkszs.^2)) |
---|
| 193 | error('Dimensions of G do not match G_blkszs.'); |
---|
| 194 | end; |
---|
| 195 | |
---|
| 196 | % mineigF is the smallest eigenvalue of F_0 |
---|
| 197 | mineigF = 0.0; |
---|
| 198 | k=0; for n=F_blkszs, |
---|
| 199 | mineigF = min(mineigF, min(eig(reshape(F(k+[1:n*n],1),n,n)))); |
---|
| 200 | k=k+n*n; % k = sum n_i*n_i |
---|
| 201 | end; |
---|
| 202 | % mineigG is the smallest eigenvalue of G_0 |
---|
| 203 | mineigG = 0.0; |
---|
| 204 | k=0; for n=G_blkszs, |
---|
| 205 | mineigG = min(mineigG, min(eig(reshape(G(k+[1:n*n],1),n,n)))); |
---|
| 206 | k=k+n*n; % k = sum n_i*n_i |
---|
| 207 | end; |
---|
| 208 | |
---|
| 209 | % eyeF is the identity |
---|
| 210 | eyeF = zeros(size(F,1),1); |
---|
| 211 | k=0; for n=F_blkszs, |
---|
| 212 | eyeF(k+[1:n*n]) = reshape(eye(n),n*n,1); % identity |
---|
| 213 | k=k+n*n; % k = sum n_i*n_i |
---|
| 214 | end; |
---|
| 215 | % eyeG is the identity |
---|
| 216 | eyeG = zeros(size(G,1),1); |
---|
| 217 | k=0; for n=G_blkszs, |
---|
| 218 | eyeG(k+[1:n*n]) = reshape(eye(n),n*n,1); % identity |
---|
| 219 | k=k+n*n; % k = sum n_i*n_i |
---|
| 220 | end; |
---|
| 221 | |
---|
| 222 | % initial x0 |
---|
| 223 | x0 = [zeros(m,1); max(-1.1*min(mineigF,mineigG), 1)]; |
---|
| 224 | |
---|
| 225 | % linear objective |
---|
| 226 | c = [zeros(m,1); 1]; |
---|
| 227 | |
---|
| 228 | % call maxdet |
---|
| 229 | [x,Z,W,ul,hist,infostr]=maxdet([F,eyeF; G,eyeG],... |
---|
| 230 | [F_blkszs(:)',G_blkszs(:)'],... |
---|
| 231 | [1 zeros(1,m+1)],1,c,x0,... |
---|
| 232 | zeros(size(F,1)+size(G,1),1),0,... |
---|
| 233 | abstol,reltol,gam,NTiters); |
---|
| 234 | |
---|
| 235 | % prepare output |
---|
| 236 | x = x(1:m); |
---|
| 237 | W = Z(size(F,1)+1:size(F,1)+size(G,1)); |
---|
| 238 | Z = Z(1:size(F,1)); |
---|
| 239 | %if ul(1)<0 |
---|
| 240 | % infostr = 'feasible'; |
---|
| 241 | %elseif ul(2)>=0 |
---|
| 242 | % infostr = 'infeasible'; |
---|
| 243 | %else |
---|
| 244 | % infostr = 'feasibility cannot be determined'; |
---|
| 245 | %end |
---|