diff --git a/Contents.m b/Contents.m index deda01f..50b2f32 100644 --- a/Contents.m +++ b/Contents.m @@ -43,4 +43,5 @@ % eyeK - Identity with respect to symmetric cones % cellK - Transforms x-solution into a cell array % -% See also conversion. \ No newline at end of file +% See also conversion. + diff --git a/blkchol2.c b/blkchol2.c index 7b5a174..c1822c5 100644 --- a/blkchol2.c +++ b/blkchol2.c @@ -113,16 +113,21 @@ void cholonBlk(double *x, double *d, mwIndex m, const mwIndex ncols, const mwInd ------------------------------------------------------- */ xkk = x[inz]; if(xkk > lb[k]){ /* now xkk > 0 */ - if(xkk < ub){ +/* ------------------------------------------------------------ + maxabs is a wrapper for the BLAS IDAMAX Fortran function. + IDAMAX finds the first element having maximum absolute + value in an array. Only call maxabs with m>1. + ------------------------------------------------------------ */ + if ((m>1) && (xkk < ub)){ ubk = maxabs(x+inz+1,m-1) / maxu; if(xkk < ubk){ /* ------------------------------------------------------------ If we need to add on diagonal, store this in (skipIr, lb(k)). ------------------------------------------------------------ */ skipIr[nskip++] = first + k; - lb[k] = ubk - xkk; /* amount added on diagonal */ - xkk = ubk; - } + lb[k] = ubk - xkk; /* amount added on diagonal */ + xkk = ubk; + } } /* -------------------------------------------------------------- Set dk = xkk, lkk = 1 (for LDL'). @@ -170,7 +175,7 @@ void cholonBlk(double *x, double *d, mwIndex m, const mwIndex ncols, const mwInd INCLUDING THE DIAGONAL ENTRY. Lir - Lir[0:nnz-1] ARE THE ROW INDICES OF THE NONZEROS OF THE FIRST COLUMN OF THE SUPERNODE. - OUTPUT PARAMETERS - + OUTPUT PARAMETERS - irInv - On return, irInv[Lir[0:nnz-1]] = nnz:-1:1, so that Lir[nnz-irInv[i]] == i The position of subscript "xij" is thus @@ -314,7 +319,7 @@ void spadd(const mwIndex *xjjc, double *xnz, const mwIndex mj, const mwIndex nj, } } -/* ************************************************************ +/* ************************************************************ PROCEDURE precorrect - Apply corrections from affecting supernode (skipping subnodes with non-positive diagonal) on supernodal diagonal block in L-factor. @@ -520,7 +525,7 @@ mwIndex blkLDL(const mwIndex neqns, const mwIndex nsuper, const mwIndex *xsuper, length[k],xsuper[k],xsuper[k+1], relind,fwsiz,fwork)) == (mwIndex)-1 ) return (mwIndex)-1; /* fwsiz too small */ - } + } /* ------------------------------------------------------------ DO DENSE CHOLESKY on the current supernode ------------------------------------------------------------ */ diff --git a/bwblkslv.c b/bwblkslv.c index 55a1a8e..f837e6d 100644 --- a/bwblkslv.c +++ b/bwblkslv.c @@ -186,7 +186,7 @@ void mexFunction(const int nlhs, mxArray *plhs[], mwIndex m,n, j, k, nsuper, inz; double *y, *fwork; const double *permPr, *b, *xsuperPr; - const mwIndex *yjc, *yir, *bjc, *bir; + const mwIndex *yjc=NULL, *yir=NULL, *bjc=NULL, *bir=NULL; mwIndex *perm, *xsuper, *iwork, *snode; jcir L; char bissparse; diff --git a/dpr1fact.c b/dpr1fact.c index 45778aa..0c92d61 100644 --- a/dpr1fact.c +++ b/dpr1fact.c @@ -365,6 +365,7 @@ char dodpr1fact(double *beta, mwIndex *perm, double *d, double t, const double * ------------------------------------------------------------ */ else{ psqrdep = 0.0; + j = 0; for(i = 0; dep[i] < m; i++) if(psqr[dep[i]] > psqrdep){ j = i; diff --git a/examples/README.md b/examples/README.md index abd62e8..1e17b6d 100644 --- a/examples/README.md +++ b/examples/README.md @@ -11,7 +11,7 @@ http://plato.asu.edu/sub/testcases.html. > Applied Mechanics and Engineering 94(1992):113-129. > DOI: [10.1016/0045-7825(92)90159-H](https://doi.org/10.1016/0045-7825(92)90159-H). - Optimal value is `5.66517e-01`. + Optimal value is `-5.66517e-01`. - **control07.mat**: Again from SDPLIB, a problem from control and system theory @@ -22,7 +22,7 @@ http://plato.asu.edu/sub/testcases.html. > Department of Mathematical and Computing Sciences, Tokyo Institute of > Technology, September, 1997. - Optimal value is `2.06251e+01`. + Optimal value is `-2.06251e+01`. - **nb.mat**: This problem is by Robert Vanderbei, it was used at the @@ -39,8 +39,12 @@ http://plato.asu.edu/sub/testcases.html. > "The reduced density matrix method for electronic structure calculations > and the role of three-index representability", October, 2003. > DOI: [10.1063/1.1636721](https://doi.org/10.1063/1.1636721). - - at http://mf.c.titech.ac.jp/mituhiro/software.html. + + and http://mf.c.titech.ac.jp/mituhiro/software.html. The reported energy is + `-75.1014` substracted by the repulsion energy `4.3656986614` and negated + due to the problem formulation. + + Optimal value is `7.946708e+01`. - **trto3.mat**: A problem by Kocvara, from single-load truss topology design. diff --git a/examples/test_sedumi.m b/examples/test_sedumi.m index 9b82b13..6226ae6 100644 --- a/examples/test_sedumi.m +++ b/examples/test_sedumi.m @@ -1,49 +1,64 @@ function test_sedumi (do_rebuild, do_exit) % TEST_SEDUMI Verify minimal functionality of SeDuMi using an example SDP. - -if (nargin < 2) - do_rebuild = false; - do_exit = false; -end +% +% Input: +% do_rebuild - Rebuild SeDuMi +% do_exit - Exit Matlab/GNU Octave on failure (suiteable for continuous +% integration CI). try old_dir = cd ('..'); - if (do_rebuild) + if ((nargin == 1) && do_rebuild) install_sedumi -rebuild else install_sedumi end cd (old_dir) - + example_path = mfilename ('fullpath'); example_path = example_path(1:max ([1, strfind(example_path, filesep ())]) - 1); - - test_problems = {'arch0.mat', 'control07.mat', 'nb.mat', ... - 'OH_2Pi_STO-6GN9r12g1T2.mat', 'trto3.mat'}; - - for i = 1:length (test_problems) + + % Name and optimal function value + test_problems = { ... + 'arch0.mat', -5.665170e-01; ... + 'control07.mat', -2.062510e+01; ... + 'nb.mat', -5.070309e-02; ... + 'OH_2Pi_STO-6GN9r12g1T2.mat', 7.946708e+01; ... + 'trto3.mat', -1.279999e+04}; + + tol = 1e-6; + + for i = 1:size (test_problems, 1) line_sep = repmat ('-', 1, 75); fprintf ('\n%s\n-- Test ''%s'' (%2d/%2d)\n%s\n\n', line_sep, ... - test_problems{i}, i, length (test_problems), line_sep); - - load (fullfile (example_path, test_problems{i}), 'At', 'b', 'c', 'K'); - [~,~,info] = sedumi (At, b, c, K); - if (info.pinf ~= 0 || info.dinf ~= 0 || info.numerr == 2) - fprintf ('\n\nSeDuMi TEST FAILED: Could not solve ''%s'' example.\n', ... - test_problems{i}); - if (do_exit) + test_problems{i,1}, i, size (test_problems, 1), line_sep); + + load (fullfile (example_path, test_problems{i,1}), 'At', 'b', 'c', 'K'); + [x, y, info] = sedumi (At, b, c, K); + cx = c' * x; + by = b' * y; + cx_relerr = abs (cx - test_problems{i,2}) / abs (test_problems{i,2}); + by_relerr = abs (by - test_problems{i,2}) / abs (test_problems{i,2}); + if ((cx_relerr >= tol) || (by_relerr >= tol) ... + || (info.pinf ~= 0) || (info.dinf ~= 0) || (info.numerr == 2)) + fprintf ('\n\tTEST FAILED\n'); + fprintf ('\n\t\tExpected = %e ( tol = %.1e)', test_problems{i,2}, tol); + fprintf ('\n\t\t c*x = %e (relerr = %.1e)', cx, cx_relerr); + fprintf ('\n\t\t b*y = %e (relerr = %.1e)\n', by, by_relerr); + if ((nargin == 2) && do_exit) exit (1); else return; end end if (info.numerr == 1) - warning ('SeDuMi reported reduced accuracy.') + fprintf ('\n\tWARNING: SeDuMi reported reduced accuracy.\n'); end + fprintf ('\n\tTEST PASSED\n'); end catch - disp ('SeDuMi TEST FAILED: Compilation or runtime error.'); - if (do_exit) + fprintf ('\n\tTEST FAILED: Compilation or runtime error.\n'); + if ((nargin == 2) && do_exit) exit (1); else return; @@ -51,3 +66,4 @@ function test_sedumi (do_rebuild, do_exit) end end + diff --git a/extractA.c b/extractA.c index ee7f683..23792c6 100644 --- a/extractA.c +++ b/extractA.c @@ -98,7 +98,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { jcir At, Apart; - mwIndex i, m, njc, ifirst, n, ynnz, blk0,blk1; + mwIndex i, m, njc, ifirst=0, n=0, ynnz, blk0, blk1; mwIndex *Ajc; const double *blkstartPr, *AjcPr; bool isblk0negative; diff --git a/fwblkslv.c b/fwblkslv.c index b449376..931cc2a 100644 --- a/fwblkslv.c +++ b/fwblkslv.c @@ -197,7 +197,7 @@ void mexFunction(const int nlhs, mxArray *plhs[], mwIndex m,n, j, k, nsuper, inz; double *y,*fwork; const double *permPr, *b, *xsuperPr; - const mwIndex *yjc, *yir, *bjc, *bir; + const mwIndex *yjc=NULL, *yir=NULL, *bjc=NULL, *bir=NULL; mwIndex *perm, *invperm, *snode, *xsuper, *iwork; jcir L; char bissparse; diff --git a/maxstep.m b/maxstep.m index fbb4663..e4ffe83 100644 --- a/maxstep.m +++ b/maxstep.m @@ -50,7 +50,7 @@ reltr = x(ix(1):ix(2)-1).*dx(ix(1):ix(2)-1)... - ddot(x(ix(2):ix(3)-1),dx,K.qblkstart); norm2 = reltr.^2 - tdet(dx,K).*auxx.tdet; - if norm2 > 0 + if all(norm2 > 0) norm2 = sqrt(norm2); end mindxq = min( (reltr - norm2)./auxx.tdet); @@ -64,4 +64,4 @@ mindxs = minpsdeig(reldx, K); mindx = min(mindx, mindxs); end -tp = 1 / max(-mindx, 1E-16); \ No newline at end of file +tp = 1 / max(-mindx, 1E-16); diff --git a/psdframeit.c b/psdframeit.c index d5a1286..f9cc2ae 100644 --- a/psdframeit.c +++ b/psdframeit.c @@ -107,9 +107,9 @@ void psdframeit(double *x, const double *frms, const double *lab, void mexFunction(const int nlhs, mxArray *plhs[], const int nrhs, const mxArray *prhs[]) { - mwIndex i,lendiag, lenfull, lenud,qsize; + mwIndex i, lendiag, lenud, qsize; double *x, *fwork; - const double *lab,*frms; + const double *lab, *frms; mwIndex *sdpNL; coneK cK; /* ------------------------------------------------------------ @@ -126,7 +126,6 @@ void mexFunction(const int nlhs, mxArray *plhs[], ------------------------------------------------------------ */ lenud = cK.rDim + cK.hDim; qsize = lenud + cK.hLen; - lenfull = cK.lpN + cK.qDim + lenud; lendiag = cK.lpN + 2 * cK.lorN + cK.rLen + cK.hLen; /* ------------------------------------------------------------ Get inputs lab,frms diff --git a/sedumi.m b/sedumi.m index f225414..5c75fa9 100644 --- a/sedumi.m +++ b/sedumi.m @@ -264,7 +264,11 @@ if N*m<100000 %Test if Ax=b is feasible at all %turn off the rank deficient warning for now - s = warning('off','MATLAB:singularMatrix'); + if (exist ('OCTAVE_VERSION', 'builtin') == 5) + s = warning('off','Octave:singular-matrix'); + else + s = warning('off','MATLAB:singularMatrix'); + endif y=[A;b']\[zeros(N,1);1]; if abs(y'*b-1) < 1e-10 && norm(A*y) < 1e-10 %Infeasibility certificate found @@ -778,8 +782,6 @@ info.err(1)=norm(x'*(origcoeff.At)-(origcoeff.b)',2)/(1+normb); %Let us get rid of the K.f part, since the free variables don't make %any difference in the cone infeasibility. - %origcoeff.K.f=0; - if origcoeff.K.f= jlast){ /* move to new z-column */ @@ -334,7 +337,8 @@ void spcpxdxd(double *z, const mwIndex *zir, const mwIndex znnz, mwIndex inz, i, icol, j, jfirst, jlast, m, nsqr, imgfirst; double *dx, *dxpi, *fworkpi; double zi; - const double *dj, *djpi, *djx, *djxpi; + const double *dj=NULL, *djpi=NULL, *djx=NULL, *djxpi=NULL; + jfirst = 0; nsqr = SQR(n); /* ------------------------------------------------------------ Partition 4*n^2 WORKING array fwork into [fwork(2*n^2), dxRows(2*n^2)]. diff --git a/sqrtinv.c b/sqrtinv.c index e74acd1..29645c4 100644 --- a/sqrtinv.c +++ b/sqrtinv.c @@ -95,7 +95,7 @@ void prpiqdivv(double *y,double *ypi, const double *q,const double *qpi, void mexFunction(const int nlhs, mxArray *plhs[], const int nrhs, const mxArray *prhs[]) { - mwIndex i,k, nk, nksqr, lenud, lendiag, diagskip; + mwIndex k, nk, nksqr, lenud, lendiag, diagskip; double *y; const double *q,*v; coneK cK; diff --git a/symfctmex.c b/symfctmex.c index 53296e5..b34e14b 100644 --- a/symfctmex.c +++ b/symfctmex.c @@ -136,7 +136,6 @@ void mexFunction(const int nlhs, mxArray *plhs[], *invp, *colcnt; mxArray *L_FIELD; const char *LFieldnames[] = {"L", "perm", "xsuper"}; - mwIndex *mwXjc, *mwXir, *mwLjc, *mwLir; /* ------------------------------------------------------------ Check for proper number of arguments ------------------------------------------------------------ */ diff --git a/trydif.m b/trydif.m index 50abd0a..0714174 100644 --- a/trydif.m +++ b/trydif.m @@ -52,7 +52,7 @@ halfxz = (x(ix(1):ix(2)-1).*z(ix(1):ix(2)-1)... + ddot(x(ix(2):ix(3)-1),z,K.qblkstart)) / 2; tmp = halfxz.^2 - detxz; - if tmp > 0 + if all(tmp > 0) lab2q = halfxz + sqrt(tmp); else lab2q = halfxz; @@ -69,4 +69,4 @@ t = 0; w = wIN; wr = wrIN; -end \ No newline at end of file +end diff --git a/urotorder.c b/urotorder.c index 01fdf6a..86b5b74 100644 --- a/urotorder.c +++ b/urotorder.c @@ -90,6 +90,7 @@ void rotorder(mwIndex *perm, double *u, mwIndex *gjc, twodouble *g, double *d, ------------------------------------------------------------ */ for(j = 0; j < n; j++) perm[j] = j; + pivk = 0; inz = 0; d[0] = 0.0; h = 1.0; for(k = 0, rowuk = u; k < n-1; k++, rowuk++){ @@ -208,6 +209,7 @@ void prpirotorder(mwIndex *perm, double *u,double *upi, mwIndex *gjc, ------------------------------------------------------------ */ for(j = 0; j < n; j++) perm[j] = j; + pivk = 0; inz = 0; d[0] = 0.0; h = 1.0; for(k = 0, rowuk = u, rowukpi = upi; k < n-1; k++, rowuk++, rowukpi++){ @@ -311,8 +313,8 @@ void mexFunction(const int nlhs, mxArray *plhs[], const int nrhs, const mxArray *prhs[]) { mxArray *myplhs[NPAROUT]; - mwIndex i,j,k, nk, nksqr, lenud, sdplen, gnnz, inz, maxKs,maxKssqr, rgnnz, hgnnz; - const double *uOld, *permOld; + mwIndex i, k, nk, nksqr, lenud, sdplen, gnnz, inz, maxKs, maxKssqr, rgnnz, hgnnz; + const double *uOld, *permOld=NULL; double *u, *d, *gjcPr, *permPr, *fwork, *fworkpi; mwIndex *perm, *gjc; double *g, *gk; diff --git a/vecsym.c b/vecsym.c index 3f116d1..991aba9 100644 --- a/vecsym.c +++ b/vecsym.c @@ -139,10 +139,8 @@ void vecsymPSD(double *y, const double *x,const mwIndex rsdpN,const mwIndex sdpN void mexFunction(const int nlhs, mxArray *plhs[], const int nrhs, const mxArray *prhs[]) { - mxArray *output_array[1], *Xk; - coneK cK; - mwIndex k, nk, nksqr, lqDim,lenfull; + mwIndex lqDim, lenfull; const double *x; double *y; diff --git a/widelen.m b/widelen.m index 2dacf29..0d49bf2 100644 --- a/widelen.m +++ b/widelen.m @@ -91,7 +91,7 @@ halfxz = (xM(ix(1):ix(2)-1).*zM(ix(1):ix(2)-1)... + ddot(xM(ix(2):ix(3)-1),zM,K.qblkstart)) / 2; tmp = halfxz.^2 - detxz; - if tmp > 0 + if all(tmp > 0) lab2q = halfxz + sqrt(tmp); else lab2q = halfxz; @@ -124,4 +124,4 @@ wr.alpha=alphaM; wr.delta=deltaM; end -wr.desc = 1; % always descent direction. \ No newline at end of file +wr.desc = 1; % always descent direction.