More improvements to the eigenvalue computations. Modified the handli…

…ng of the EIGS call in maxeigK and minpsdeig based on feedback from Jon Dattoro; now SeDuMi handles convergence failures by reverting to the standard eigenvalue computation.

eigK.m

@@ -50,89 +50,77 @@
 % 02110-1301, USA
 %
 
-% The existence of rsdpN is code for 'is this the internal format?'
-is_int = isfield(K,'rsdpN');
-if is_int,
-    N = K.l+2*length(K.q)+sum(K.s);
-else    
+if isfield(K,'rsdpN'),
+    % The internal SeDuMi cone format.
+    is_int = true;
+    nf = 0;
+    nl = K.l;
+    nq = length(K.q);
+    nr = 0;
+    ns = length(K.s);
+    nrsdp = K.rsdpN;
+    N = nl+2*length(K.q)+sum(K.s);
+else
+    % The external SeDuMi cone format.
+    is_int = false;
     N = 0;
-    if isfield(K,'l'), N = N + K.l; end
-    if isfield(K,'q'), N = N + 2*length(K.q); end
-    if isfield(K,'r'), N = N + 2*length(K.r); end
-    if isfield(K,'s'), N = N + sum(K.s); end
+    if isfield(K,'f'), nf = K.f; else nf = 0; end
+    if isfield(K,'l'), nl = K.l; N = N + nl; else nl = 0; end
+    if isfield(K,'q'), nq = length(K.q); N = N + 2 * nq; else nq = 0; end
+    if isfield(K,'r'), nr = length(K.r); N = N + 2 * nr; else nr = 0; end
+    if isfield(K,'s'), ns = length(K.s); N = N + sum(K.s); else ns = 0; end
 end
 li = 0;
-xi = 0;
+xi = nf;
 lab = zeros(N,1);
-if ~is_int && isfield(K,'f'),
-    xi = xi + K.f;
-end
-if isfield(K,'l') && K.l > 0,
-    li(li+1:li+K.l) = x(xi+1:xi+K.l);
-    xi = xi+K.l;
-    li = li+K.l;
-end
-if isfield(K,'q') && ~isempty(K.q),
+li(li+1:li+nl) = x(xi+1:xi+nl);
+xi = xi + nl;
+li = li + nl;
+if nq,
     tmp = sqrt(0.5);
     if is_int,
         % Internal version: all of the x0 values are at the front, and the
         % vectors are stacked in order after that.
         zi = xi;
-        xi = xi + length(K.q);
-        for k = 1:length(K.q)
-            kk = K.q(k) - 1;
-            x0 = x(zi+k);
+        xi = xi + nq;
+        for i = 1:nq,
+            kk = K.q(i) - 1;
+            x0 = x(zi+i);
             lab(li+1:li+2) = tmp*(x0+[-1;+1]*norm(x(xi+1:xi+kk)));
             xi = xi + kk;
             li = li + 2;
         end
     else
-        for k = 1:length(K.q)
-            kk = K.q(k);
+        for i = 1:length(K.q)
+            kk = K.q(i);
             x0 = x(xi+1);
             lab(li+1:li+2) = tmp*(x0+[-1;+1]*norm(x(xi+2:xi+kk)));
             xi = xi + kk;
             li = li + 2;
         end
     end
 end
-if ~is_int && isfield(K,'r') && ~isempty(K.r),
-    % Rotated cones are not encountered internally, so only the standard
-    % ordering is implemented, and a simpler formula than the one found
-    % in eigK.c is being used, since cancellation error is not a concern.
-    tmp = 0.5;
-    for k = 1:length(K.r)
-        kk = K.r(k);
-        x1 = xx(xi+1);
-        x2 = xx(xi+2);
-        lab(li+1:li+2) = tmp*(x1+x2+[-1;+1]*norm([x1-x2;2*x(xi+3:xi+kk)]));
-        xi = xi + kk;
-        li = li + 2;
-    end
+for i = 1:nr,
+    % Only the external format need be implemented here
+    ki = K.r(i);
+    x1 = xx(xi+1);
+    x2 = xx(xi+2);
+    lab(li+1:li+2) = 0.5*(x1+x2+[-1;+1]*norm([x1-x2;2*x(xi+3:xi+ki)]));
+    xi = xi + ki;
+    li = li + 2;
 end
-if isfield(K,'s') && ~isempty(K.s),
-    Ks = K.s;
-    Kq = K.s .* K.s;
-    nc = length(Ks);
-    % When used internally, Hermitian terms are broken apart into real and
-    % imaginary halves.
-    if is_int,
-        nr = K.rsdpN;
-    else
-        nr = nc;
-    end
-    for i = 1 : nc,
-        ki = Ks(i);
-        qi = Kq(i);
-        XX = x(xi+1:xi+qi); 
+for i = 1 : ns,
+    ki = K.s(i);
+    qi = ki * ki;
+    XX = x(xi+1:xi+qi); 
+    xi = xi + qi;
+    if i > nrsdp,
+        XX = XX + 1i*x(xi+1:xi+qi); 
         xi = xi + qi;
-        if i > nr,
-            XX = XX + 1i*x(xi+1:xi+qi); 
-            xi = xi + qi;
-        end
-        XX = reshape(XX,ki,ki);
-        lab(li+1:li+ki) = 0.5 * eig( XX + XX' );
-        li = li + ki;
     end
+    XX = reshape(XX,ki,ki);
+    XX = XX + XX';
+    lab(li+1:li+ki) = 0.5*eig( XX );
+    li = li + ki;
 end
 

maxeigK.m

@@ -1,6 +1,6 @@
 function lab = maxeigK(x,K)
 
-% lab = mineigK(x,K)
+% lab = maxeigK(x,K)
 %
 % MAXEIGK  Computes the maximum eigenvalue of a vector x with respect to a 
 %          self-dual homogenous cone K.
@@ -9,59 +9,95 @@
 
 % New function by Michael C. Grant
 % Copyright (C) 2013 Michael C. Grant.
+%
+% This program is free software; you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation; either version 2 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program; if not, write to the Free Software
+% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
+% 02110-1301, USA
+%
 
-xi = 0;
-if isfield(K,'f'),
-    xi = xi + K.f;
-end
-if isfield(K,'l') && K.l > 0,
-    lab = max(x(xi+1:xi+K.l));
-    xi = xi+K.l;
+% The existence of rsdpN is code for 'is this the internal format?'
+if isfield(K,'rsdpN'),
+    % The internal SeDuMi cone format.
+    is_int = true;
+    nf = 0;
+    nl = K.l;
+    nq = length(K.q);
+    nr = 0;
+    ns = length(K.s);
+    nrsdp = K.rsdpN;
 else
-    lab = -Inf;
+    % The external SeDuMi cone format.
+    is_int = false;
+    if isfield(K,'f'), nf = K.f; else nf = 0; end
+    if isfield(K,'l'), nl = K.l; else nl = 0; end
+    if isfield(K,'q'), nq = length(K.q); else nq = 0; end
+    if isfield(K,'r'), nr = length(K.r); else nr = 0; end
+    if isfield(K,'s'), ns = length(K.s); else ns = 0; end
 end
-if isfield(K,'q') && ~isempty(K.q),
-    scl = sqrt(0.5);
-    for k = 1:length(K.q)
-        kk = K.q(k);
-        lab = max(lab,scl*(x(xi+1)+norm(x(xi+2:xi+kk))));
-        xi = xi + kk;
+xi = nf;
+lab = max([-Inf;x(xi+1:xi+nl)]);
+xi = xi + nl;
+if nq,
+    tmp = sqrt(0.5);
+    if is_int,
+        % Internal version: all of the x0 values are at the front, and the
+        % vectors are stacked in order after that.
+        zi = xi;
+        xi = xi + nq;
+        for i = 1:nq,
+            kk = K.q(i) - 1;
+            x0 = x(zi+i);
+            lab = max(lab,tmp*(x0+norm(x(xi+1:xi+kk))));
+            xi = xi + kk;
+        end
+    else
+        for i = 1:length(K.q)
+            kk = K.q(i);
+            x0 = x(xi+1);
+            lab = max(lab,tmp*(x0+norm(x(xi+2:xi+kk))));
+            xi = xi + kk;
+        end
     end
 end
-if isfield(K,'r') && ~isempty(K.r),
-    for k = 1:length(K.r)
-        kk = K.r(k);
-        x1 = xx(xi+1);
-        x2 = xx(xi+2);
-        lab = max(lab,0.5*(x1+x2+norm([x1-x2;2*x(xi+3:xi+kk)])));
-    end
+for i = 1:nr,
+    % Only the external format need be implemented here
+    ki = K.r(i);
+    x1 = xx(xi+1);
+    x2 = xx(xi+2);
+    lab = max(lab,0.5*(x1+x2+norm([x1-x2;2*x(xi+3:xi+ki)])));
+    xi = xi + ki;
 end
-if isfield(K,'s') && ~isempty(K.s),
-    Ks = K.s;
-    Kq = K.s .* K.s;
-    nc = length(Ks);
-    OPTS.disp=0;
-    % When used internally, Hermitian terms are broken apart into real and
-    % imaginary halves, so we need to catch this.
-    if isfield(K,'rsdpN'),
-        nr = K.rsdpN;
-    else
-        nr = nc;
-    end
-    for i = 1 : nc,
-        ki = Ks(i);
-        qi = Kq(i);
-        XX = x(xi+1:xi+qi); xi=xi+qi;
-        if i > nr,
-            XX = XX + 1i*x(xi+1:xi+qi); xi=xi+qi;
+if ns,
+    for i = 1 : ns,
+        ki = K.s(i);
+        qi = ki * ki;
+        XX = x(xi+1:xi+qi); 
+        xi = xi + qi;
+        if i > nrsdp,
+            XX = XX + 1i*x(xi+1:xi+qi); 
+            xi = xi + qi;
         end
         XX = reshape(XX,ki,ki);
         XX = XX + XX';
-        if ki > 500
-            lab=max(lab,0.5*eigs(XX,1,'LA',OPTS));
+        if ki > 500,
+            if nnz(XX) < 0.1 * numel(XX), XX = sparse(XX); end
+            [v,val,flag] = eigs(XX,1,'LA',struct('issym',true)); %#ok
+            if flag, val = max(eig(XX)); end
         else
-            lab=max(lab,0.5*max(eig(XX)));
+            val = max(eig(XX));
         end
+        lab = max(lab,0.5*val);
     end
 end
 

minpsdeig.m

@@ -52,11 +52,8 @@
 Kq = Ks .* Ks;
 nr = K.rsdpN;
 nc = length(Ks);
-N  = sum(Kq) + sum(Kq(nr+1:end));
-ux = zeros(N,1);
-xi = length(x) - N;
+xi = length(x) - sum(Kq) - sum(Kq(nr+1:end));
 eigv = zeros(nc,1);
-OPTS.disp=0;
 for i = 1 : nc,
     ki = Ks(i);
     qi = Kq(i);
@@ -68,8 +65,10 @@
     end
     XX = reshape(XX,ki,ki);
     XX = XX + XX';
-    if ki > 500
-        eigv(i) = eigs(XX,1,'SA',OPTS);
+    if ki > 500,
+        if nnz(XX) < 0.1 * numel(XX), XX = sparse(XX); end
+        [v,eigv(i),flag] = eigs(XX,1,'SA',struct('issym',true)); %#ok
+        if flag, eigv(i) = min(eig(XX)); end
     else
         eigv(i) = min(eig(XX));
     end

sedumi.m

@@ -176,14 +176,6 @@
 %    measures as defined in the Seventh DIMACS Challenge. For more details
 %    see the User Guide.
 %
-%    (13) pars.free   By default, pars.free=1, and all free variables are
-%                     collected into a single second-order cone block. If
-%                     pars.free=0, then they are each split into the
-%                     difference off two nonnegative variables. If
-%                     pars.free=2, then SeDuMi uses the split approach if
-%                     the rest of the problem is an LP, and the
-%                     second-order cone approach otherwise.
-%
 % Bug reports can be submitted at http://sedumi.mcmaster.ca.
 %
 % See also mat, vec, cellK, eyeK, eigK
@@ -548,7 +540,7 @@
             break
         end
     elseif (by > 0) && (abs(1+feasratio) < 0.05) && (R.b0*y0 < 0.5)
-        if maxeigK(qreshape(Amul(A,dense,y,1),1,K),K) <= pars.eps * by
+        if maxeigK(Amul(A,dense,y,1),K) <= pars.eps * by
             STOP = 3;                   % Means Farkas solution found !
             break
         end
@@ -624,16 +616,15 @@
 % Determine infeasibility
 % ------------------------------------------------------------
 pinf = norm(x0*b-Ax);
-z = qreshape(Ay-x0*c,1,K);
-dinf = maxeigK(z,K);
+dinf = maxeigK(Ay-x0*c,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 = maxeigK(qreshape(Ay,1,K),K);
+        ddirinf = maxeigK(Ay,K);
         if cx < 0.0
             reldirinf = pdirinf / (-cx);
         else