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 |
---|