53 #pragma GCC diagnostic ignored "-Wstrict-prototypes"
54 #include "sdpa_call.h"
55 #pragma GCC diagnostic warning "-Wstrict-prototypes"
57 #include "blockmemshell/memory.h"
59 #include "scip/pub_misc.h"
63 #define GAPTOLCHANGE 1
64 #define FEASTOLCHANGE 1
65 #define PENALTYBOUNDTOL 1E-3
68 #define INFEASFEASTOLCHANGE 0.1
69 #define INFEASMINFEASTOL 1E-9
71 #define MIN_LAMBDASTAR 1e0
72 #define MAX_LAMBDASTAR 1e8
73 #define LAMBDASTAR_FACTOR 1e0
74 #define LAMBDASTAR_TWOPOINTS TRUE
75 #define LAMBDASTAR_THRESHOLD 1e1
76 #define LAMBDASTAR_LOW 1.5
77 #define LAMBDASTAR_HIGH 1e5
79 #define MIN_PENALTYPARAM 1e5
80 #define MAX_PENALTYPARAM 1e12
81 #define PENALTYPARAM_FACTOR 1e1
82 #define MAX_MAXPENALTYPARAM 1e15
83 #define MAXPENALTYPARAM_FACTOR 1e6
86 #define BMS_CALL(x) do \
90 SCIPerrorMessage("No memory in function call.\n"); \
91 return SCIP_NOMEMORY; \
97 #define CHECK_IF_SOLVED(sdpisolver) do \
99 if (!(sdpisolver->solved)) \
101 SCIPerrorMessage("Tried to access solution information for SDP %d ahead of solving!\n", sdpisolver->sdpcounter); \
102 return SCIP_LPERROR; \
108 #define CHECK_IF_SOLVED_BOOL(sdpisolver) do \
110 if (!(sdpisolver->solved)) \
112 SCIPerrorMessage("Tried to access solution information for SDP %d ahead of solving!\n", sdpisolver->sdpcounter); \
120 struct SCIP_SDPiSolver
122 SCIP_MESSAGEHDLR* messagehdlr;
128 int* inputtosdpamapper;
131 int* sdpatoinputmapper;
132 SCIP_Real* fixedvarsval;
133 SCIP_Real fixedvarsobjcontr;
145 SCIP_Real sdpsolverfeastol;
151 SCIP_Real lambdastar;
168 assert( sdpisolver != NULL );
169 assert( lb < ub + sdpisolver->feastol );
171 return (ub-lb <= sdpisolver->epsilon);
174 #define isFixed(sdpisolver,lb,ub) (ub-lb <= sdpisolver->epsilon)
182 SCIP_Real penaltyparam,
190 int* sdpconstnblocknonz,
194 SCIP_Real** sdpconstval,
196 int** sdpnblockvarnonz,
206 int* blockindchanges,
220 char phase_string[15];
223 assert( feastol != NULL );
227 SCIP_Real* solvector;
229 SCIP_Bool infeasible;
232 BMS_CALL( BMSallocBufferMemoryArray(sdpisolver->bufmem, &solvector, nvars) );
233 nvarspointer = nvars;
235 assert( nvarspointer == nvars );
238 SCIP_CALL(
SCIPsdpSolcheckerCheck(sdpisolver->bufmem, nvars, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
239 sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval,
240 indchanges, nremovedinds, blockindchanges, nlpcons, noldlpcons, lplhs, lprhs, rownactivevars, lpnnonz, lprow, lpcol, lpval,
241 solvector, sdpisolver->feastol, sdpisolver->epsilon, &infeasible) );
243 BMSfreeBufferMemoryArray(sdpisolver->bufmem, &solvector);
247 SCIPdebugMessage(
"Solution feasible for SDPA but outside feasibility tolerance, changing SDPA feasibility tolerance from %f to %f\n",
254 sdpisolver->sdpa->setParameterEpsilonDash(*feastol);
256 #ifdef SCIP_MORE_DEBUG
257 sdpisolver->sdpa->printParameters(stdout);
259 sdpisolver->sdpa->setInitPoint(
false);
260 #ifdef SDPA_RESETPARAMS
261 sdpisolver->sdpa->resetParameters();
263 sdpisolver->sdpa->initializeSolve();
265 sdpisolver->sdpa->solve();
268 sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
269 sdpisolver->nsdpcalls += 1;
273 sdpisolver->sdpa->getPhaseString((
char*)phase_string);
274 SCIPdebugMessage(
"SDPA solving finished with status %s (primal and dual here are switched in contrast to our formulation)\n", phase_string);
279 sdpisolver->solved = FALSE;
280 SCIPmessagePrintInfo(sdpisolver->messagehdlr,
"SDPA failed to reach required feasibility tolerance! \n");
310 return "Primal-dual Interior Point Solver for SDPs developed by K. Fujisawa et al. (sdpa.sourceforge.net)";
323 assert( sdpisolver != NULL );
324 return (
void*) sdpisolver->sdpa;
356 SCIP_MESSAGEHDLR* messagehdlr,
361 assert( sdpisolver != NULL );
362 assert( blkmem != NULL );
363 assert( bufmem != NULL );
365 SCIPdebugMessage(
"Calling SCIPsdpiCreate \n");
367 BMS_CALL( BMSallocBlockMemory(blkmem, sdpisolver) );
369 (*sdpisolver)->messagehdlr = messagehdlr;
370 (*sdpisolver)->blkmem = blkmem;
371 (*sdpisolver)->bufmem = bufmem;
374 (*sdpisolver)->sdpa = NULL;
376 (*sdpisolver)->nvars = 0;
377 (*sdpisolver)->nactivevars = 0;
378 (*sdpisolver)->inputtosdpamapper = NULL;
379 (*sdpisolver)->sdpatoinputmapper = NULL;
380 (*sdpisolver)->fixedvarsval = NULL;
381 (*sdpisolver)->fixedvarsobjcontr = 0.0;
382 (*sdpisolver)->nvarbounds = 0;
383 (*sdpisolver)->varboundpos = NULL;
384 (*sdpisolver)->solved = FALSE;
385 (*sdpisolver)->timelimit = FALSE;
386 (*sdpisolver)->sdpcounter = 0;
387 (*sdpisolver)->niterations = 0;
388 (*sdpisolver)->nsdpcalls = 0;
390 (*sdpisolver)->epsilon = 1e-9;
391 (*sdpisolver)->gaptol = 1e-4;
392 (*sdpisolver)->feastol = 1e-6;
393 (*sdpisolver)->sdpsolverfeastol = 1e-6;
395 (*sdpisolver)->sdpinfo = FALSE;
406 assert( sdpisolver != NULL );
407 assert( *sdpisolver != NULL );
409 SCIPdebugMessage(
"Freeing SDPISolver\n");
412 if ( (*sdpisolver)->sdpa != NULL)
413 delete (*sdpisolver)->sdpa;
415 BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->varboundpos, 2 * (*sdpisolver)->nactivevars);
417 BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->inputtosdpamapper, (*sdpisolver)->nvars);
419 BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->sdpatoinputmapper, (*sdpisolver)->nactivevars);
421 BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->fixedvarsval, (*sdpisolver)->nvars - (*sdpisolver)->nactivevars);
423 BMSfreeBlockMemory((*sdpisolver)->blkmem, sdpisolver);
433 assert( sdpisolver != NULL );
435 sdpisolver->sdpcounter++;
445 assert( sdpisolver != NULL );
447 SCIPdebugMessage(
"Resetting counter of SDP-Interface from %d to 0.\n", sdpisolver->sdpcounter);
448 sdpisolver->sdpcounter = 0;
485 int* sdpconstnblocknonz,
489 SCIP_Real** sdpconstval,
491 int** sdpnblockvarnonz,
501 int* blockindchanges,
507 int* lprownactivevars,
518 return SCIPsdpiSolverLoadAndSolveWithPenalty(sdpisolver, 0.0, TRUE, FALSE, nvars, obj, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
519 sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
520 nremovedinds, blockindchanges, nremovedblocks, nlpcons, noldlpcons, lplhs, lprhs, lprownactivevars, lpnnonz, lprow, lpcol, lpval, start,
521 startsettings, timelimit, NULL, NULL);
545 SCIP_Real penaltyparam,
556 int* sdpconstnblocknonz,
560 SCIP_Real** sdpconstval,
562 int** sdpnblockvarnonz,
572 int* blockindchanges,
589 SCIP_Bool* penaltybound
594 SCIP_Real* sdpavarbounds;
607 #ifdef SCIP_MORE_DEBUG
612 char phase_string[15];
615 assert( sdpisolver != NULL );
616 assert( penaltyparam > -1 * sdpisolver->epsilon );
617 assert( penaltyparam < sdpisolver->epsilon || ( feasorig != NULL ) );
619 assert( obj != NULL );
620 assert( lb != NULL );
621 assert( ub != NULL );
622 assert( nsdpblocks >= 0 );
623 assert( nsdpblocks == 0 || sdpblocksizes != NULL );
624 assert( nsdpblocks == 0 || sdpnblockvars != NULL );
625 assert( sdpconstnnonz >= 0 );
626 assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstnblocknonz != NULL );
627 assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstrow != NULL );
628 assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstcol != NULL );
629 assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstval != NULL );
630 assert( sdpnnonz >= 0 );
631 assert( nsdpblocks == 0 || sdpnblockvarnonz != NULL );
632 assert( nsdpblocks == 0 || sdpvar != NULL );
633 assert( nsdpblocks == 0 || sdprow != NULL );
634 assert( nsdpblocks == 0 || sdpcol != NULL );
635 assert( nsdpblocks == 0 || sdpval != NULL );
636 assert( nsdpblocks == 0 || indchanges != NULL );
637 assert( nsdpblocks == 0 || nremovedinds != NULL );
638 assert( nsdpblocks == 0 || blockindchanges != NULL );
639 assert( 0 <= nremovedblocks && nremovedblocks <= nsdpblocks );
640 assert( nlpcons >= 0 );
641 assert( noldlpcons >= nlpcons );
642 assert( nlpcons == 0 || lplhs != NULL );
643 assert( nlpcons == 0 || lprhs != NULL );
644 assert( nlpcons == 0 || rownactivevars != NULL );
645 assert( lpnnonz >= 0 );
646 assert( nlpcons == 0 || lprow != NULL );
647 assert( nlpcons == 0 || lpcol != NULL );
648 assert( nlpcons == 0 || lpval != NULL );
650 sdpisolver->niterations = 0;
651 sdpisolver->nsdpcalls = 0;
654 if ( timelimit <= 0.0 )
656 sdpisolver->solved = FALSE;
657 sdpisolver->timelimit = TRUE;
661 sdpisolver->timelimit = FALSE;
673 if ( penaltyparam < sdpisolver->epsilon )
674 SCIPdebugMessage(
"Inserting Data into SDPA for SDP (%d) \n", ++sdpisolver->sdpcounter);
676 SCIPdebugMessage(
"Inserting Data again into SDPA for SDP (%d) \n", sdpisolver->sdpcounter);
679 sdpisolver->penalty = (penaltyparam < sdpisolver->epsilon) ? FALSE : TRUE;
680 sdpisolver->rbound = rbound;
684 BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->inputtosdpamapper), sdpisolver->nvars, nvars) );
685 BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->sdpatoinputmapper), sdpisolver->nactivevars, nvars) );
686 BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->fixedvarsval), sdpisolver->nvars - sdpisolver->nactivevars, nvars) );
688 oldnactivevars = sdpisolver->nactivevars;
689 sdpisolver->nvars = nvars;
690 sdpisolver->nactivevars = 0;
694 sdpisolver->fixedvarsobjcontr = 0.0;
695 for (i = 0; i < nvars; i++)
697 if (
isFixed(sdpisolver, lb[i], ub[i]) )
699 sdpisolver->fixedvarsobjcontr += obj[i] * lb[i];
700 sdpisolver->fixedvarsval[nfixedvars] = lb[i];
702 sdpisolver->inputtosdpamapper[i] = -nfixedvars;
703 SCIPdebugMessage(
"Fixing variable %d locally to %f for SDP %d in SDPA\n", i, lb[i], sdpisolver->sdpcounter);
707 sdpisolver->sdpatoinputmapper[sdpisolver->nactivevars] = i;
708 sdpisolver->nactivevars++;
709 sdpisolver->inputtosdpamapper[i] = sdpisolver->nactivevars;
710 #ifdef SCIP_MORE_DEBUG
711 SCIPdebugMessage(
"Variable %d becomes variable %d for SDP %d in SDPA\n", i, sdpisolver->inputtosdpamapper[i], sdpisolver->sdpcounter);
715 assert( sdpisolver->nactivevars + nfixedvars == sdpisolver->nvars );
719 sdpisolver->fixedvarsobjcontr = 0.0;
722 BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->fixedvarsval), nvars, nfixedvars) );
723 BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->sdpatoinputmapper), nvars, sdpisolver->nactivevars) );
726 if ( sdpisolver->sdpa != 0 )
729 delete sdpisolver->sdpa;
730 sdpisolver->sdpa =
new SDPA();
733 sdpisolver->sdpa =
new SDPA();
734 assert( sdpisolver->sdpa != 0 );
739 sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_STABLE_BUT_SLOW);
743 feastol = FEASTOLCHANGE * FEASTOLCHANGE * sdpisolver->sdpsolverfeastol;
744 SCIPdebugMessage(
"Start solving process with stable settings\n");
748 sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_UNSTABLE_BUT_FAST);
749 sdpisolver->sdpa->setParameterEpsilonStar(sdpisolver->gaptol);
750 sdpisolver->sdpa->setParameterEpsilonDash(sdpisolver->sdpsolverfeastol);
751 feastol = sdpisolver->sdpsolverfeastol;
752 SCIPdebugMessage(
"Start solving process with fast settings\n");
756 sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_DEFAULT);
758 sdpisolver->sdpa->setParameterEpsilonStar(
GAPTOLCHANGE * sdpisolver->gaptol);
759 sdpisolver->sdpa->setParameterEpsilonDash(
FEASTOLCHANGE * sdpisolver->sdpsolverfeastol);
761 SCIPdebugMessage(
"Start solving process with medium settings\n");
765 SCIPdebugMessage(
"Unknown setting for start-settings: %d!\n", startsettings); \
768 sdpisolver->sdpa->setParameterLowerBound(-1e20);
773 sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
775 sdpisolver->sdpa->setParameterUpperBound(1e8);
778 sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
780 #ifdef SCIP_MORE_DEBUG
781 sdpisolver->sdpa->printParameters(stdout);
785 sdpisolver->nvarbounds = 0;
786 BMS_CALL( BMSallocBufferMemoryArray(sdpisolver->bufmem, &sdpavarbounds, 2 * sdpisolver->nactivevars) );
788 if ( sdpisolver->nactivevars != oldnactivevars )
790 if ( sdpisolver->varboundpos == NULL )
792 BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->varboundpos), 2 * sdpisolver->nactivevars) );
796 BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->varboundpos), 2 * oldnactivevars, 2 * sdpisolver->nactivevars) );
799 assert( sdpisolver->varboundpos != NULL );
801 for (i = 0; i < sdpisolver->nactivevars; i++)
803 assert( 0 <= sdpisolver->sdpatoinputmapper[i] && sdpisolver->sdpatoinputmapper[i] < nvars );
806 sdpavarbounds[sdpisolver->nvarbounds] = lb[sdpisolver->sdpatoinputmapper[i]];
807 sdpisolver->varboundpos[sdpisolver->nvarbounds] = -(i + 1);
808 (sdpisolver->nvarbounds)++;
812 sdpavarbounds[sdpisolver->nvarbounds] = ub[sdpisolver->sdpatoinputmapper[i]];
813 sdpisolver->varboundpos[sdpisolver->nvarbounds] = +(i + 1);
814 (sdpisolver->nvarbounds)++;
821 BMS_CALL( BMSallocBufferMemoryArray(sdpisolver->bufmem, &rowmapper, 2*noldlpcons) );
826 for (i = 0; i < noldlpcons; i++)
828 if ( rownactivevars[i] >= 2 )
832 rowmapper[2*i] = pos;
839 rowmapper[2*i + 1] = pos;
843 rowmapper[2*i + 1] = -1;
850 rowmapper[2*i + 1] = -1;
854 assert( nlpineqs <= 2*nlpcons );
860 if ( penaltyparam >= sdpisolver->epsilon && rbound )
861 sdpisolver->nvarbounds++;
863 if ( sdpisolver->sdpinfo )
864 sdpisolver->sdpa->setDisplay(stdout);
866 sdpisolver->sdpa->setDisplay(0);
868 #ifdef SCIP_MORE_DEBUG
869 FILE* fpOut = fopen(
"output.tmp",
"w");
872 sdpisolver->sdpa->setResultFile(fpOut);
876 if ( penaltyparam < sdpisolver->epsilon )
877 sdpisolver->sdpa->inputConstraintNumber((
long long) sdpisolver->nactivevars);
879 sdpisolver->sdpa->inputConstraintNumber((
long long) sdpisolver->nactivevars + 1);
882 sdpisolver->sdpa->inputBlockNumber((
long long) ((nlpineqs + sdpisolver->nvarbounds > 0) ?
883 nsdpblocks - nremovedblocks + 1 : nsdpblocks - nremovedblocks));
886 for (block = 0; block < nsdpblocks; block++)
888 if ( blockindchanges[block] >= 0 )
890 SCIPdebugMessage(
"adding block %d to SDPA as block %d with size %d\n",
891 block, block - blockindchanges[block] + 1, sdpblocksizes[block] - nremovedinds[block]);
892 sdpisolver->sdpa->inputBlockSize((
long long) block - blockindchanges[block] + 1,
893 (
long long) sdpblocksizes[block] - nremovedinds[block]);
894 sdpisolver->sdpa->inputBlockType((
long long) block - blockindchanges[block] + 1, SDPA::SDP);
897 if ( nlpineqs + sdpisolver->nvarbounds > 0 )
900 sdpisolver->sdpa->inputBlockSize((
long long) nsdpblocks - nremovedblocks + 1, (
long long) -(nlpineqs + sdpisolver->nvarbounds));
901 sdpisolver->sdpa->inputBlockType((
long long) nsdpblocks - nremovedblocks + 1, SDPA::LP);
902 SCIPdebugMessage(
"adding LP block to SDPA as block %d with size %d\n", nsdpblocks - nremovedblocks + 1,
903 -(nlpineqs + sdpisolver->nvarbounds));
906 sdpisolver->sdpa->initializeUpperTriangleSpace();
909 for (i = 0; i < sdpisolver->nactivevars; i++)
914 sdpisolver->sdpa->inputCVec((
long long) i + 1, obj[sdpisolver->sdpatoinputmapper[i]]);
915 #ifdef SCIP_MORE_DEBUG
916 SCIPdebugMessage(
"inserting objective %f for variable %d which became variable %d in SDPA\n", obj[sdpisolver->sdpatoinputmapper[i]],
917 sdpisolver->sdpatoinputmapper[i], i+1);
920 if ( penaltyparam >= sdpisolver->epsilon )
921 sdpisolver->sdpa->inputCVec((
long long) sdpisolver->nactivevars + 1, penaltyparam);
926 sdpisolver->sdpa->setInitPoint(
true);
928 sdpisolver->sdpa->setInitPoint(
false);
936 assert( nsdpblocks > 0 );
937 assert( sdpnblockvarnonz != NULL );
938 assert( sdpnblockvars != NULL );
939 assert( sdpcol != NULL );
940 assert( sdprow != NULL );
941 assert( sdpval != NULL );
942 assert( sdpvar != NULL );
943 assert( indchanges != NULL );
944 assert( nremovedinds != NULL );
946 for (block = 0; block < nsdpblocks; block++)
949 if ( blockindchanges[block] == -1 )
951 #ifdef SCIP_MORE_DEBUG
952 SCIPdebugMessage(
" -> building block %d, which becomes block %d in SDPA (%d)\n", block, block - blockindchanges[block] + 1,sdpisolver->sdpcounter);
955 for (blockvar = 0; blockvar < sdpnblockvars[block]; blockvar++)
957 v = sdpisolver->inputtosdpamapper[sdpvar[block][blockvar]];
959 #ifdef SCIP_MORE_DEBUG
960 SCIPdebugMessage(
" -> adding coefficient matrix for variable %d which becomes variable %d in SDPA (%d)\n",
961 sdpvar[block][blockvar], v, sdpisolver->sdpcounter);
967 for (k = 0; k < sdpnblockvarnonz[block][blockvar]; k++)
970 assert( indchanges[block][sdprow[block][blockvar][k]] > -1 && indchanges[block][sdpcol[block][blockvar][k]] > -1 );
972 assert( indchanges[block][sdprow[block][blockvar][k]] <= sdprow[block][blockvar][k]);
973 assert( indchanges[block][sdpcol[block][blockvar][k]] <= sdpcol[block][blockvar][k]);
975 assert( 0 <= sdprow[block][blockvar][k] && sdprow[block][blockvar][k] < sdpblocksizes[block] );
976 assert( 0 <= sdpcol[block][blockvar][k] && sdpcol[block][blockvar][k] < sdpblocksizes[block] );
979 #ifdef SCIP_MORE_DEBUG
980 SCIPdebugMessage(
" -> adding nonzero %g at (%d,%d) (%d)\n",
981 sdpval[block][blockvar][k],
982 sdpcol[block][blockvar][k] - indchanges[block][sdpcol[block][blockvar][k]] + 1,
983 sdprow[block][blockvar][k] - indchanges[block][sdprow[block][blockvar][k]] + 1,
984 sdpisolver->sdpcounter);
987 sdpisolver->sdpa->inputElement((
long long) v, (
long long) block - blockindchanges[block] + 1,
988 (
long long) sdpcol[block][blockvar][k] - indchanges[block][sdpcol[block][blockvar][k]] + 1,
989 (
long long) sdprow[block][blockvar][k] - indchanges[block][sdprow[block][blockvar][k]] + 1,
990 sdpval[block][blockvar][k], checkinput);
995 if ( penaltyparam >= sdpisolver->epsilon )
997 #ifdef SCIP_MORE_DEBUG
998 SCIPdebugMessage(
" -> adding coefficient matrix for penalty variable r in SDPA (%d)\n", sdpisolver->sdpcounter);
1000 for (i = 0; i < sdpblocksizes[block] - nremovedinds[block]; i++)
1002 #ifdef SCIP_MORE_DEBUG
1003 SCIPdebugMessage(
" -> adding nonzero 1.0 at (%d,%d) (%d)\n", i + 1, i + 1, sdpisolver->sdpcounter);
1006 sdpisolver->sdpa->inputElement((
long long) sdpisolver->nactivevars + 1, (
long long) block - blockindchanges[block] + 1,
1007 (
long long) i + 1, (
long long) i + 1, 1.0, checkinput);
1014 if ( sdpconstnnonz > 0 )
1016 assert( nsdpblocks > 0 );
1017 assert( sdpconstnblocknonz!= NULL );
1018 assert( sdpconstcol != NULL );
1019 assert( sdpconstrow != NULL );
1020 assert( sdpconstval != NULL );
1022 for (block = 0; block < nsdpblocks; block++)
1024 if ( blockindchanges[block] == -1 )
1026 #ifdef SCIP_MORE_DEBUG
1027 SCIPdebugMessage(
" -> building block %d (%d)\n", block + 1, sdpisolver->sdpcounter);
1029 for (k = 0; k < sdpconstnblocknonz[block]; k++)
1032 assert( indchanges[block][sdpconstrow[block][k]] > -1 && indchanges[block][sdpconstcol[block][k]] > -1 );
1034 assert( indchanges[block][sdpconstrow[block][k]] <= sdpconstrow[block][k]);
1035 assert( indchanges[block][sdpconstcol[block][k]] <= sdpconstcol[block][k]);
1037 assert (0 <= sdpconstrow[block][k] && sdpconstrow[block][k] < sdpblocksizes[block]);
1038 assert (0 <= sdpconstcol[block][k] && sdpconstcol[block][k] < sdpblocksizes[block]);
1041 #ifdef SCIP_MORE_DEBUG
1042 SCIPdebugMessage(
" -> adding constant nonzero %g at (%d,%d) (%d)\n", sdpconstval[block][k],
1043 sdpconstcol[block][k] - indchanges[block][sdpconstcol[block][k]] + 1,
1044 sdpconstrow[block][k] - indchanges[block][sdpconstrow[block][k]] + 1,
1045 sdpisolver->sdpcounter);
1047 sdpisolver->sdpa->inputElement((
long long) 0, (
long long) block - blockindchanges[block] + 1,
1048 (
long long) sdpconstcol[block][k] - indchanges[block][sdpconstcol[block][k]] + 1,
1049 (
long long) sdpconstrow[block][k] - indchanges[block][sdpconstrow[block][k]] + 1,
1050 sdpconstval[block][k], checkinput);
1055 #ifdef SCIP_MORE_DEBUG
1056 SCIPdebugMessage(
" -> building LP-block %d (%d)\n", nsdpblocks - nremovedblocks + 1, sdpisolver->sdpcounter);
1060 for (i = 0; i < lpnnonz; i++)
1062 assert( 0 <= lprow[i] && lprow[i] < noldlpcons );
1063 assert( 0 <= lpcol[i] && lpcol[i] < nvars );
1064 assert( REALABS(lpval[i]) > sdpisolver->epsilon );
1067 if ( sdpisolver->inputtosdpamapper[lpcol[i]] > 0 )
1070 assert( rownactivevars[lprow[i]] > 0 );
1071 if ( rownactivevars[lprow[i]] > 1 )
1073 if ( lprow[i] > lastrow )
1077 if ( penaltyparam >= sdpisolver->epsilon )
1080 if ( rowmapper[2*lastrow] > -1 )
1082 #ifdef SCIP_MORE_DEBUG
1083 SCIPdebugMessage(
" -> adding nonzero 1.0 at (%d,%d) for penalty variable r in SDPA (%d)\n",
1084 rowmapper[2*lastrow], rowmapper[2*lastrow], sdpisolver->sdpcounter);
1088 sdpisolver->sdpa->inputElement((
long long) sdpisolver->nactivevars + 1, (
long long) nsdpblocks - nremovedblocks + 1,
1089 (
long long) rowmapper[2*lastrow], (
long long) rowmapper[2*lastrow], 1.0, checkinput);
1093 if ( rowmapper[2*lastrow + 1] > -1 )
1095 #ifdef SCIP_MORE_DEBUG
1096 SCIPdebugMessage(
" -> adding nonzero 1.0 at (%d,%d) for penalty variable r in SDPA (%d)\n",
1097 rowmapper[2*lastrow + 1], rowmapper[2*lastrow + 1], sdpisolver->sdpcounter);
1101 sdpisolver->sdpa->inputElement((
long long) sdpisolver->nactivevars + 1, (
long long) nsdpblocks - nremovedblocks + 1,
1102 (
long long) rowmapper[2*lastrow + 1], (
long long) rowmapper[2*lastrow + 1], 1.0, checkinput);
1107 if ( rowmapper[2*lastrow] > -1 )
1109 #ifdef SCIP_MORE_DEBUG
1110 SCIPdebugMessage(
" -> adding nonzero %g at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1111 lpval[i], rowmapper[2*lastrow], rowmapper[2*lastrow], lpcol[i], sdpisolver->inputtosdpamapper[lpcol[i]], sdpisolver->sdpcounter);
1114 sdpisolver->sdpa->inputElement((
long long) sdpisolver->inputtosdpamapper[lpcol[i]], (
long long) nsdpblocks - nremovedblocks + 1,
1115 (
long long) rowmapper[2*lastrow], (
long long) rowmapper[2*lastrow], lpval[i], checkinput);
1118 if ( rowmapper[2*lastrow + 1] > -1 )
1120 #ifdef SCIP_MORE_DEBUG
1121 SCIPdebugMessage(
" -> adding nonzero %g at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1122 -1 * lpval[i], rowmapper[2*lastrow + 1], rowmapper[2*lastrow + 1], lpcol[i], sdpisolver->inputtosdpamapper[lpcol[i]], sdpisolver->sdpcounter);
1126 sdpisolver->sdpa->inputElement((
long long) sdpisolver->inputtosdpamapper[lpcol[i]], (
long long) nsdpblocks - nremovedblocks + 1,
1127 (
long long) rowmapper[2*lastrow + 1], (
long long) rowmapper[2*lastrow + 1], -1 * lpval[i], checkinput);
1135 for (i = 0; i < nlpcons; i++)
1140 #ifdef SCIP_MORE_DEBUG
1141 SCIPdebugMessage(
" -> adding lhs %g at (%d,%d) (%d)\n", lplhs[i], lpconsind, lpconsind, sdpisolver->sdpcounter);
1144 sdpisolver->sdpa->inputElement((
long long) 0, (
long long) nsdpblocks - nremovedblocks + 1, (
long long) lpconsind,
1145 (
long long) lpconsind, lplhs[i], checkinput);
1152 #ifdef SCIP_MORE_DEBUG
1153 SCIPdebugMessage(
" -> adding lhs (originally rhs) %g at (%d,%d) (%d)\n", -1 * lprhs[i], lpconsind, lpconsind, sdpisolver->sdpcounter);
1157 sdpisolver->sdpa->inputElement((
long long) 0, (
long long) nsdpblocks - nremovedblocks + 1, (
long long) lpconsind,
1158 (
long long) lpconsind, -1 * lprhs[i], checkinput);
1163 assert( lpconsind == nlpineqs + 1 );
1166 #ifdef SCIP_MORE_DEBUG
1169 for (i = 0; i < lpnnonz; i++)
1172 if ( sdpisolver->inputtosdpamapper[lpcol[i]] > 0 )
1174 if ( rownactivevars[lprow[i]] > 1 )
1176 if ( lprow[i] > lastrow )
1181 SCIPmessagePrintInfo(sdpisolver->messagehdlr,
" <= %f\n", lprhs[ind]);
1183 SCIPmessagePrintInfo(sdpisolver->messagehdlr,
"\n");
1188 SCIPmessagePrintInfo(sdpisolver->messagehdlr,
"%f <= ", lplhs[ind]);
1190 SCIPmessagePrintInfo(sdpisolver->messagehdlr,
"+ %f <x%d> ", lpval[i], lpcol[i]);
1197 SCIPmessagePrintInfo(sdpisolver->messagehdlr,
" <= %f\n", lprhs[ind]);
1199 SCIPmessagePrintInfo(sdpisolver->messagehdlr,
"\n");
1201 assert( ind == nlpcons - 1 );
1206 BMSfreeBufferMemoryArray(sdpisolver->bufmem, &rowmapper);
1210 for (i = 0; i < ((penaltyparam < sdpisolver->epsilon) || (! rbound) ? sdpisolver->nvarbounds : sdpisolver->nvarbounds - 1); i++)
1212 assert( 0 < abs(sdpisolver->varboundpos[i]) && abs(sdpisolver->varboundpos[i] <= sdpisolver->nactivevars) );
1215 if ( sdpisolver->varboundpos[i] < 0 )
1219 sdpisolver->sdpa->inputElement((
long long) -sdpisolver->varboundpos[i], (
long long) nsdpblocks - nremovedblocks + 1,
1220 (
long long) nlpineqs + 1 + i, (
long long) nlpineqs + 1 + i, 1.0, checkinput);
1222 if ( REALABS(sdpavarbounds[i]) > sdpisolver->epsilon )
1225 #ifdef SCIP_MORE_DEBUG
1226 SCIPdebugMessage(
" -> adding lower bound %g at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1227 sdpavarbounds[i], nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpatoinputmapper[-sdpisolver->varboundpos[i] - 1],
1228 -sdpisolver->varboundpos[i], sdpisolver->sdpcounter);
1230 sdpisolver->sdpa->inputElement((
long long) 0, (
long long) nsdpblocks - nremovedblocks + 1, (
long long) nlpineqs + 1 + i,
1231 (
long long) nlpineqs + 1 + i, sdpavarbounds[i], checkinput);
1236 #ifdef SCIP_MORE_DEBUG
1237 SCIPdebugMessage(
" -> adding lower bound 0 at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1238 nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpatoinputmapper[-sdpisolver->varboundpos[i] - 1],
1239 -sdpisolver->varboundpos[i], sdpisolver->sdpcounter);
1249 sdpisolver->sdpa->inputElement((
long long) sdpisolver->varboundpos[i], (
long long) nsdpblocks - nremovedblocks + 1,
1250 (
long long) nlpineqs + 1 + i, (
long long) nlpineqs + 1 + i, -1.0, checkinput);
1252 if ( REALABS(sdpavarbounds[i]) > sdpisolver->epsilon )
1255 #ifdef SCIP_MORE_DEBUG
1256 SCIPdebugMessage(
" -> adding upper bound %g at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1257 sdpavarbounds[i], nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpatoinputmapper[sdpisolver->varboundpos[i] - 1],
1258 sdpisolver->varboundpos[i], sdpisolver->sdpcounter);
1260 sdpisolver->sdpa->inputElement((
long long) 0, (
long long) nsdpblocks - nremovedblocks + 1, (
long long) nlpineqs + 1 + i,
1261 (
long long) nlpineqs + 1 + i, -sdpavarbounds[i], checkinput);
1266 #ifdef SCIP_MORE_DEBUG
1267 SCIPdebugMessage(
" -> adding upper bound 0 at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1268 0, nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpatoinputmapper[sdpisolver->varboundpos[i] - 1],
1269 sdpisolver->varboundpos[i]);
1275 if ( penaltyparam >= sdpisolver->epsilon && rbound )
1278 sdpisolver->sdpa->inputElement((
long long) sdpisolver->nactivevars + 1, (
long long) nsdpblocks - nremovedblocks + 1,
1279 (
long long) nlpineqs + 1 + i, (
long long) nlpineqs + 1 + i, 1.0, checkinput);
1280 #ifdef SCIP_MORE_DEBUG
1281 SCIPdebugMessage(
" -> adding lower bound r >= 0 at (%d,%d) in SDPA (%d)\n", nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpcounter);
1286 BMSfreeBufferMemoryArray(sdpisolver->bufmem, &sdpavarbounds);
1289 sdpisolver->sdpa->initializeUpperTriangle();
1290 sdpisolver->sdpa->initializeSolve();
1293 if ( start != NULL && penaltyparam < sdpisolver->epsilon )
1295 SCIPdebugMessage(
"Starting with a previous solution is not yet tested for the interface, only x-vector is given, not y and Z");
1296 for (i = 1; i <= sdpisolver->nactivevars; i++)
1297 sdpisolver->sdpa->inputInitXVec((
long long) i, start[sdpisolver->sdpatoinputmapper[i] - 1]);
1299 else if ( penaltyparam >= sdpisolver->epsilon )
1300 SCIPdebugMessage(
"Skipping insertion of starting point, as this is not yet supported for penalty formulation.\n");
1302 #ifdef SCIP_DEBUG_PRINTTOFILE
1304 sdpisolver->sdpa->writeInputSparse(const_cast<char*>(
"sdpa.dat-s"), const_cast<char*>(
"%+8.3e"));
1305 sdpisolver->sdpa->writeInitSparse(const_cast<char*>(
"sdpa.ini-s"), const_cast<char*>(
"%+8.3e"));
1308 SCIPdebugMessage(
"Calling SDPA solve (SDP: %d)\n", sdpisolver->sdpcounter);
1309 sdpisolver->sdpa->solve();
1310 sdpisolver->solved = TRUE;
1313 sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1314 sdpisolver->nsdpcalls += 1;
1318 sdpisolver->sdpa->getPhaseString((
char*)phase_string);
1319 SCIPdebugMessage(
"SDPA solving finished with status %s (primal and dual here are switched in contrast to our formulation)\n", phase_string);
1325 else if ( penaltyparam >= sdpisolver->epsilon )
1330 SCIP_CALL(
checkFeastolAndResolve(sdpisolver, penaltyparam, nvars, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
1331 sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
1332 nremovedinds, blockindchanges, nlpcons, noldlpcons, lplhs, lprhs, rownactivevars, lpnnonz, lprow, lpcol, lpval, &feastol) );
1339 SCIPdebugMessage(
"Numerical troubles -- solving SDP %d again ...\n", sdpisolver->sdpcounter);
1342 sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_DEFAULT);
1343 sdpisolver->sdpa->setParameterEpsilonStar(
GAPTOLCHANGE * sdpisolver->gaptol);
1344 sdpisolver->sdpa->setParameterEpsilonDash(
FEASTOLCHANGE * sdpisolver->sdpsolverfeastol);
1345 sdpisolver->sdpa->setParameterLowerBound(-1e20);
1348 sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
1350 sdpisolver->sdpa->setParameterUpperBound(1e8);
1353 sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
1355 #ifdef SCIP_MORE_DEBUG
1356 sdpisolver->sdpa->printParameters(stdout);
1358 sdpisolver->sdpa->setInitPoint(
false);
1359 #ifdef SDPA_RESETPARAMS
1360 sdpisolver->sdpa->resetParameters();
1362 sdpisolver->sdpa->initializeSolve();
1364 sdpisolver->sdpa->solve();
1365 sdpisolver->solved = TRUE;
1368 sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1369 sdpisolver->nsdpcalls += 1;
1377 sdpisolver->sdpa->getPhaseString((
char*)phase_string);
1378 SCIPdebugMessage(
"SDPA solving finished with status %s (primal and dual here are switched in contrast to our formulation)\n", phase_string);
1383 SCIP_CALL(
checkFeastolAndResolve(sdpisolver, penaltyparam, nvars, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
1384 sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
1385 nremovedinds, blockindchanges, nlpcons, noldlpcons, lplhs, lprhs, rownactivevars, lpnnonz, lprow, lpcol, lpval, &feastol) );
1392 SCIPdebugMessage(
"Numerical troubles -- solving SDP %d again^2 ...\n", sdpisolver->sdpcounter);
1395 sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_STABLE_BUT_SLOW);
1398 sdpisolver->sdpa->setParameterLowerBound(-1e20);
1401 sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
1403 sdpisolver->sdpa->setParameterUpperBound(1e8);
1406 sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
1408 #ifdef SCIP_MORE_DEBUG
1409 sdpisolver->sdpa->printParameters(stdout);
1411 sdpisolver->sdpa->setInitPoint(
false);
1412 #ifdef SDPA_RESETPARAMS
1413 sdpisolver->sdpa->resetParameters();
1415 sdpisolver->sdpa->initializeSolve();
1417 sdpisolver->sdpa->solve();
1418 sdpisolver->solved = TRUE;
1421 sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1422 sdpisolver->nsdpcalls += 1;
1430 sdpisolver->sdpa->getPhaseString((
char*)phase_string);
1431 SCIPdebugMessage(
"SDPA solving finished with status %s (primal and dual here are switched in constrast to our formulation)\n", phase_string);
1436 SCIP_CALL(
checkFeastolAndResolve(sdpisolver, penaltyparam, nvars, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
1437 sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
1438 nremovedinds, blockindchanges, nlpcons, noldlpcons, lplhs, lprhs, rownactivevars, lpnnonz, lprow, lpcol, lpval, &feastol) );
1441 #ifdef SCIP_MORE_DEBUG
1442 (void) fclose(fpOut);
1446 if ( penaltyparam >= sdpisolver->epsilon )
1453 SCIP_Real trace = 0.0;
1456 assert( (sdpisolver->nactivevars + 1 == sdpisolver->sdpa->getConstraintNumber()) );
1458 sdpasol = sdpisolver->sdpa->getResultXVec();
1461 *feasorig = (sdpasol[sdpisolver->nactivevars] < sdpisolver->feastol);
1465 if ( ! *feasorig && penaltybound != NULL )
1469 SCIPdebugMessage(
"Solution not feasible in original problem, r = %f\n", sdpasol[sdpisolver->nactivevars]);
1474 nblockssdpa = (int) sdpisolver->sdpa->getBlockNumber();
1475 for (b = 1; b <= nblockssdpa; b++)
1478 X = sdpisolver->sdpa->getResultYMat((
long long) b);
1479 nrow = (int) sdpisolver->sdpa->getBlockSize((
long long) b);
1480 assert( nrow >= 0 );
1483 if ( sdpisolver->sdpa->getBlockType((
long long) b) == SDPA::LP )
1486 for (i = 0; i < nrow - sdpisolver->nvarbounds; i++)
1492 for (i = 0; i < nrow; i++)
1493 trace += X[i + i*nrow];
1500 if ( penaltybound != NULL )
1501 *penaltybound = TRUE;
1502 SCIPdebugMessage(
"Tr(X) = %f == %f = Gamma, penalty formulation not exact, Gamma should be increased or problem is infeasible\n",
1503 trace, penaltyparam);
1505 else if ( penaltybound != NULL )
1506 *penaltybound = FALSE;
1532 assert( sdpisolver != NULL );
1533 return sdpisolver->solved;
1546 SDPA::PhaseType phasetype;
1548 assert( sdpisolver != NULL );
1549 assert( sdpisolver->sdpa != NULL);
1552 phasetype = sdpisolver->sdpa->getPhaseValue();
1554 if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1563 SCIP_Bool* primalfeasible,
1564 SCIP_Bool* dualfeasible
1567 SDPA::PhaseType phasetype;
1569 assert( sdpisolver != NULL );
1570 assert( sdpisolver->sdpa != NULL );
1571 assert( primalfeasible != NULL );
1572 assert( dualfeasible != NULL );
1575 phasetype = sdpisolver->sdpa->getPhaseValue();
1577 switch ( phasetype )
1580 *primalfeasible = TRUE;
1581 *dualfeasible = TRUE;
1584 *primalfeasible = TRUE;
1585 *dualfeasible = TRUE;
1587 case SDPA::pFEAS_dINF:
1588 *primalfeasible = TRUE;
1589 *dualfeasible = FALSE;
1591 case SDPA::pINF_dFEAS:
1592 *primalfeasible = FALSE;
1593 *dualfeasible = TRUE;
1596 *primalfeasible = TRUE;
1597 *dualfeasible = FALSE;
1598 SCIPdebugMessage(
"SDPA stopped because dual objective became smaller than lower bound\n");
1601 *primalfeasible = FALSE;
1602 *dualfeasible = TRUE;
1603 SCIPdebugMessage(
"SDPA stopped because primal objective became bigger than upper bound\n");
1606 SCIPerrorMessage(
"SDPA doesn't know if primal and dual solutions are feasible\n");
1607 return SCIP_LPERROR;
1619 SDPA::PhaseType phasetype;
1621 assert( sdpisolver != NULL );
1622 assert( sdpisolver->sdpa != NULL );
1625 phasetype = sdpisolver->sdpa->getPhaseValue();
1627 if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::pdINF )
1629 SCIPdebugMessage(
"SDPA doesn't know if primal problem is unbounded");
1632 else if ( phasetype == SDPA::pFEAS_dINF )
1634 else if ( phasetype == SDPA::pUNBD )
1636 SCIPdebugMessage(
"SDPA was stopped because primal objective became bigger than upper bound");
1649 SDPA::PhaseType phasetype;
1651 assert( sdpisolver != NULL );
1652 assert( sdpisolver->sdpa != NULL );
1655 phasetype = sdpisolver->sdpa->getPhaseValue();
1657 if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1659 SCIPdebugMessage(
"SDPA doesn't know if primal problem is infeasible");
1662 else if ( phasetype == SDPA::pINF_dFEAS )
1664 else if ( phasetype == SDPA::dUNBD )
1666 SCIPdebugMessage(
"SDPA was stopped because dual objective became smaller than lower bound");
1679 SDPA::PhaseType phasetype;
1681 assert( sdpisolver != NULL );
1682 assert( sdpisolver->sdpa != NULL );
1685 phasetype = sdpisolver->sdpa->getPhaseValue();
1687 if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1689 SCIPdebugMessage(
"SDPA doesn't know if primal problem is feasible");
1692 else if ( phasetype == SDPA::pFEAS_dINF || phasetype == SDPA::pdOPT || phasetype == SDPA::pFEAS || phasetype == SDPA::pdFEAS )
1694 else if ( phasetype == SDPA::dUNBD )
1696 SCIPdebugMessage(
"SDPA was stopped because dual objective became smaller than lower bound");
1709 SDPA::PhaseType phasetype;
1711 assert( sdpisolver != NULL );
1712 assert( sdpisolver->sdpa != NULL);
1715 phasetype = sdpisolver->sdpa->getPhaseValue();
1717 if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1719 SCIPdebugMessage(
"SDPA doesn't know if dual problem is unbounded");
1722 else if ( phasetype == SDPA::pINF_dFEAS )
1724 else if ( phasetype == SDPA::dUNBD )
1726 SCIPdebugMessage(
"SDPA was stopped because dual objective became smaller than lower bound");
1739 SDPA::PhaseType phasetype;
1741 assert( sdpisolver != NULL );
1742 assert( sdpisolver->sdpa != NULL);
1745 phasetype = sdpisolver->sdpa->getPhaseValue();
1747 if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::pdINF )
1749 SCIPdebugMessage(
"SDPA doesn't know if dual problem is infeasible");
1752 else if ( phasetype == SDPA::pFEAS_dINF )
1754 else if ( phasetype == SDPA::pUNBD )
1756 SCIPdebugMessage(
"SDPA was stopped because primal objective became bigger than upper bound");
1769 SDPA::PhaseType phasetype;
1771 assert( sdpisolver != NULL );
1772 assert( sdpisolver->sdpa != NULL);
1775 phasetype = sdpisolver->sdpa->getPhaseValue();
1777 if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1779 SCIPdebugMessage(
"SDPA doesn't know if primal problem is feasible");
1782 else if ( phasetype == SDPA::pINF_dFEAS || phasetype == SDPA::pdOPT || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
1784 else if ( phasetype == SDPA::dUNBD )
1786 SCIPdebugMessage(
"SDPA was stopped because dual objective became smaller than lower bound");
1798 SDPA::PhaseType phasetype;
1800 assert( sdpisolver != NULL );
1801 assert( sdpisolver->sdpa != NULL );
1804 phasetype = sdpisolver->sdpa->getPhaseValue();
1806 if ( phasetype == SDPA::pdOPT )
1817 SDPA::PhaseType phasetype;
1819 assert( sdpisolver != NULL );
1820 assert( sdpisolver->sdpa != NULL );
1823 phasetype = sdpisolver->sdpa->getPhaseValue();
1825 if ( phasetype == SDPA::pUNBD )
1836 SDPA::PhaseType phasetype;
1838 assert( sdpisolver != NULL );
1839 assert( sdpisolver->sdpa != NULL);
1842 phasetype = sdpisolver->sdpa->getPhaseValue();
1844 if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
1846 if ( sdpisolver->sdpa->getParameterMaxIteration() == sdpisolver->sdpa->getIteration() )
1858 assert( sdpisolver != NULL );
1860 return sdpisolver->timelimit;
1877 SDPA::PhaseType phasetype;
1879 assert( sdpisolver != NULL );
1880 assert( sdpisolver->sdpa != NULL );
1882 if ( sdpisolver->sdpa == NULL || (! sdpisolver->solved) )
1885 phasetype = sdpisolver->sdpa->getPhaseValue();
1887 if ( phasetype == SDPA::pdOPT || phasetype == SDPA::pFEAS_dINF || phasetype == SDPA::pINF_dFEAS )
1889 if ( phasetype == SDPA::pdINF )
1891 if ( phasetype == SDPA::pUNBD)
1893 if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
1904 SDPA::PhaseType phasetype;
1906 assert( sdpisolver != NULL );
1907 assert( sdpisolver->sdpa != NULL );
1909 if ( ! sdpisolver->solved )
1912 phasetype = sdpisolver->sdpa->getPhaseValue();
1914 if ( phasetype == SDPA::pdOPT )
1926 SDPA::PhaseType phasetype;
1928 assert( sdpisolver != NULL );
1929 assert( sdpisolver->sdpa != NULL );
1931 if ( sdpisolver->timelimit )
1934 if ( ! sdpisolver->solved )
1937 phasetype = sdpisolver->sdpa->getPhaseValue();
1941 if (
SCIPsdpiSolverIsConverged(sdpisolver) || phasetype == SDPA::pUNBD || phasetype == SDPA::pINF_dFEAS || phasetype == SDPA::pFEAS_dINF )
1953 SCIPdebugMessage(
"Not implemented yet\n");
1955 return SCIP_LPERROR;
1964 assert( sdpisolver != NULL );
1965 assert( sdpisolver->sdpa != NULL );
1966 assert( objval != NULL );
1969 *objval = sdpisolver->sdpa->getPrimalObj();
1972 SCIP_Real primalval = sdpisolver->sdpa->getDualObj();
1973 SCIP_Real gap = (REALABS(*objval - primalval) / (0.5 * (REALABS(primalval) + REALABS(*objval))));
1974 if ( gap > sdpisolver->gaptol )
1975 SCIPdebugMessage(
"Attention: got objective value (before adding values of fixed variables) of %f in SCIPsdpiSolverGetObjval, "
1976 "but primal objective is %f with duality gap %f!\n", *objval, primalval, gap );
1980 *objval += sdpisolver->fixedvarsobjcontr;
2000 assert( sdpisolver != NULL );
2001 assert( sdpisolver->sdpa != NULL );
2002 assert( dualsollength != NULL );
2005 if ( objval != NULL )
2007 *objval = sdpisolver->sdpa->getPrimalObj();
2010 SCIP_Real primalval = sdpisolver->sdpa->getDualObj();
2011 SCIP_Real gap = (REALABS(*objval - primalval) / (0.5 * (REALABS(primalval) + REALABS(*objval))));
2012 if ( gap > sdpisolver->gaptol )
2014 SCIPdebugMessage(
"Attention: got objective value (before adding values of fixed variables) of %f in SCIPsdpiSolverGetSol, "
2015 "but primal objective is %f with duality gap %f!\n", *objval, primalval, gap );
2020 *objval += sdpisolver->fixedvarsobjcontr;
2023 if ( *dualsollength > 0 )
2025 assert( dualsol != NULL );
2026 if ( *dualsollength < sdpisolver->nvars )
2028 SCIPdebugMessage(
"The given array in SCIPsdpiSolverGetSol only had length %d, but %d was needed", *dualsollength, sdpisolver->nvars);
2029 *dualsollength = sdpisolver->nvars;
2035 assert( (sdpisolver->penalty && sdpisolver->nactivevars + 1 == sdpisolver->sdpa->getConstraintNumber()) ||
2036 sdpisolver->nactivevars == sdpisolver->sdpa->getConstraintNumber() );
2037 sdpasol = sdpisolver->sdpa->getResultXVec();
2039 for (v = 0; v < sdpisolver->nvars; v++)
2041 if ( sdpisolver->inputtosdpamapper[v] > 0 )
2044 dualsol[v] = sdpasol[sdpisolver->inputtosdpamapper[v] - 1];
2049 assert( -sdpisolver->inputtosdpamapper[v] <= sdpisolver->nvars - sdpisolver->nactivevars );
2050 dualsol[v] = sdpisolver->fixedvarsval[(-1 * sdpisolver->inputtosdpamapper[v]) - 1];
2077 assert( sdpisolver != NULL );
2078 assert( sdpisolver->sdpa != NULL );
2079 assert( lbvars != NULL );
2080 assert( ubvars != NULL );
2081 assert( arraylength != NULL );
2082 assert( *arraylength >= 0 );
2086 if ( *arraylength < sdpisolver->nvars )
2088 *arraylength = sdpisolver->nvars;
2089 SCIPdebugMessage(
"Insufficient length of array in SCIPsdpiSolverGetPrimalBoundVars (gave %d, needed %d)\n", *arraylength, sdpisolver->nvars);
2094 for (i = 0; i < sdpisolver->nvars; i++)
2101 if ( sdpisolver->nvarbounds == 0 )
2103 SCIPdebugMessage(
"Asked for PrimalBoundVars, but there were no variable bounds in sdpa, returning zero vector !");
2108 lpblockind = (int) sdpisolver->sdpa->getBlockNumber();
2109 assert( sdpisolver->sdpa->getBlockType((
long long) lpblockind) == SDPA::LP );
2110 nlpcons = (int) sdpisolver->sdpa->getBlockSize((
long long) lpblockind);
2111 assert( nlpcons >= 0 );
2113 X = sdpisolver->sdpa->getResultYMat((
long long) lpblockind);
2116 assert( sdpisolver->nvarbounds <= 2 * sdpisolver->nvars || (sdpisolver->nvarbounds <= 2 * sdpisolver->nvars + 1 && sdpisolver->penalty ) );
2118 for (i = 0; i < ((sdpisolver->penalty) ? sdpisolver->nvarbounds - 1 : sdpisolver->nvarbounds); i++)
2120 if ( sdpisolver->varboundpos[i] < 0 )
2124 lbvars[sdpisolver->sdpatoinputmapper[-1 * sdpisolver->varboundpos[i] -1]] = X[nlpcons - sdpisolver->nvarbounds + i];
2130 ubvars[sdpisolver->sdpatoinputmapper[+1 * sdpisolver->varboundpos[i] - 1]] = X[nlpcons - sdpisolver->nvarbounds + i];
2143 assert( sdpisolver != NULL );
2144 assert( sdpisolver->sdpa != NULL );
2145 assert( iterations != NULL );
2147 *iterations = sdpisolver->niterations;
2158 assert( sdpisolver != NULL );
2159 assert( sdpisolver->sdpa != NULL );
2160 assert( calls != NULL );
2162 *calls = sdpisolver->nsdpcalls;
2173 assert( sdpisolver != NULL );
2174 assert( usedsetting != NULL );
2177 *usedsetting = sdpisolver->usedsetting;
2218 assert( sdpisolver != NULL );
2219 assert( dval != NULL );
2224 *dval = sdpisolver->epsilon;
2227 *dval = sdpisolver->gaptol;
2230 *dval = sdpisolver->feastol;
2233 *dval = sdpisolver->sdpsolverfeastol;
2237 SCIPdebugMessage(
"Parameter SCIP_SDPPAR_PENALTYPARAM not used by SDPA");
2240 *dval = sdpisolver->objlimit;
2243 *dval = sdpisolver->lambdastar;
2246 return SCIP_PARAMETERUNKNOWN;
2259 assert( sdpisolver != NULL );
2264 sdpisolver->epsilon = dval;
2265 SCIPdebugMessage(
"Setting sdpisolver epsilon to %f.\n", dval);
2268 sdpisolver->gaptol = dval;
2269 SCIPdebugMessage(
"Setting sdpisolver gaptol to %f.\n", dval);
2272 sdpisolver->feastol = dval;
2273 SCIPdebugMessage(
"Setting sdpisolver feastol to %f.\n", dval);
2276 sdpisolver->sdpsolverfeastol = dval;
2277 SCIPdebugMessage(
"Setting sdpisolver sdpsolverfeastol to %f.\n", dval);
2280 SCIPdebugMessage(
"Parameter SCIP_SDPPAR_PENALTYPARAM not used by SDPA");
2283 SCIPdebugMessage(
"Setting sdpisolver objlimit to %f.\n", dval);
2284 sdpisolver->objlimit = dval;
2287 SCIPdebugMessage(
"Setting sdpisolver lambdastar parameter to %f.\n", dval);
2288 sdpisolver->lambdastar = dval;
2291 return SCIP_PARAMETERUNKNOWN;
2304 assert( sdpisolver != NULL );
2309 *ival = (int) sdpisolver->sdpinfo;
2310 SCIPdebugMessage(
"Getting sdpisolver information output (%d).\n", *ival);
2313 return SCIP_PARAMETERUNKNOWN;
2326 assert( sdpisolver != NULL );
2331 sdpisolver->sdpinfo = (SCIP_Bool) ival;
2332 SCIPdebugMessage(
"Setting sdpisolver information output (%d).\n", ival);
2335 return SCIP_PARAMETERUNKNOWN;
2349 assert( sdpisolver != NULL );
2354 #ifdef LAMBDASTAR_TWOPOINTS
2375 sdpisolver->lambdastar = compval;
2376 SCIPdebugMessage(
"Setting lambdastar to %f.\n", compval);
2386 SCIP_Real* penaltyparam
2391 assert( sdpisolver != NULL );
2392 assert( penaltyparam != NULL );
2408 SCIPdebugMessage(
"Setting penaltyparameter to %f.\n", compval);
2409 *penaltyparam = compval;
2417 SCIP_Real penaltyparam,
2418 SCIP_Real* maxpenaltyparam
2423 assert( sdpisolver != NULL );
2424 assert( maxpenaltyparam != NULL );
2430 *maxpenaltyparam = compval;
2431 SCIPdebugMessage(
"Setting maximum penaltyparameter to %f.\n", compval);
2460 SCIPdebugMessage(
"Not implemented yet\n");
2461 return SCIP_LPERROR;
2470 assert( fname != NULL );
2472 sdpisolver->sdpa->writeInputSparse(const_cast<char*>(fname), const_cast<char*>(
"%8.3f"));
const char * SCIPsdpiSolverGetSolverName(void)
static SCIP_RETCODE checkFeastolAndResolve(SCIP_SDPISOLVER *sdpisolver, SCIP_Real penaltyparam, int nvars, SCIP_Real *lb, SCIP_Real *ub, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nlpcons, int noldlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int *rownactivevars, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *feastol)
SCIP_RETCODE SCIPsdpiSolverSetRealpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, SCIP_Real dval)
SCIP_RETCODE SCIPsdpiSolverFree(SCIP_SDPISOLVER **sdpisolver)
SCIP_Real SCIPsdpiSolverInfinity(SCIP_SDPISOLVER *sdpisolver)
void * SCIPsdpiSolverGetSolverPointer(SCIP_SDPISOLVER *sdpisolver)
#define INFEASFEASTOLCHANGE
SCIP_RETCODE SCIPsdpiSolverIgnoreInstability(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *success)
enum SCIP_SDPSolverSetting SCIP_SDPSOLVERSETTING
SCIP_Bool SCIPsdpiSolverIsPrimalUnbounded(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetObjval(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *objval)
#define PENALTYPARAM_FACTOR
const char * SCIPsdpiSolverGetSolverDesc(void)
SCIP_Bool SCIPsdpiSolverIsInfinity(SCIP_SDPISOLVER *sdpisolver, SCIP_Real val)
interface methods for specific SDP-solvers
SCIP_RETCODE SCIPsdpiSolverResetCounter(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetRealpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, SCIP_Real *dval)
SCIP_Bool SCIPsdpiSolverIsObjlimExc(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetSol(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *objval, SCIP_Real *dualsol, int *dualsollength)
#define MAXPENALTYPARAM_FACTOR
SCIP_Bool SCIPsdpiSolverFeasibilityKnown(SCIP_SDPISOLVER *sdpisolver)
#define MAX_MAXPENALTYPARAM
SCIP_RETCODE SCIPsdpiSolverComputeLambdastar(SCIP_SDPISOLVER *sdpisolver, SCIP_Real maxguess)
SCIP_Real SCIPsdpiSolverGetDefaultSdpiSolverFeastol(void)
SCIP_Bool SCIPsdpiSolverIsPrimalFeasible(SCIP_SDPISOLVER *sdpisolver)
#define LAMBDASTAR_THRESHOLD
SCIP_RETCODE SCIPsdpiSolverSettingsUsed(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPSOLVERSETTING *usedsetting)
SCIP_RETCODE SCIPsdpiSolverGetIntpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, int *ival)
SCIP_Bool SCIPsdpiSolverWasSolved(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpSolcheckerCheck(BMS_BUFMEM *bufmem, int nvars, SCIP_Real *lb, SCIP_Real *ub, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nlpcons, int noldlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int *rownactivevars, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *solvector, SCIP_Real feastol, SCIP_Real epsilon, SCIP_Bool *infeasible)
SCIP_RETCODE SCIPsdpiSolverCreate(SCIP_SDPISOLVER **sdpisolver, SCIP_MESSAGEHDLR *messagehdlr, BMS_BLKMEM *blkmem, BMS_BUFMEM *bufmem)
#define LAMBDASTAR_FACTOR
#define CHECK_IF_SOLVED_BOOL(sdpisolver)
checks a given SDP solution for feasibility
static SCIP_Bool isFixed(SCIP_SDPISOLVER *sdpisolver, SCIP_Real lb, SCIP_Real ub)
int SCIPsdpiSolverGetDefaultSdpiSolverNpenaltyIncreases(void)
SCIP_RETCODE SCIPsdpiSolverIncreaseCounter(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetSolFeasibility(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *primalfeasible, SCIP_Bool *dualfeasible)
SCIP_RETCODE SCIPsdpiSolverWriteSDP(SCIP_SDPISOLVER *sdpisolver, const char *fname)
SCIP_Bool SCIPsdpiSolverIsAcceptable(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverIsOptimal(SCIP_SDPISOLVER *sdpisolver)
#define CHECK_IF_SOLVED(sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetSdpCalls(SCIP_SDPISOLVER *sdpisolver, int *calls)
SCIP_RETCODE SCIPsdpiSolverComputeMaxPenaltyparam(SCIP_SDPISOLVER *sdpisolver, SCIP_Real penaltyparam, SCIP_Real *maxpenaltyparam)
SCIP_RETCODE SCIPsdpiSolverGetPrimalBoundVars(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *lbvars, SCIP_Real *ubvars, int *arraylength)
SCIP_Bool SCIPsdpiSolverIsDualInfeasible(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverSetIntpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, int ival)
SCIP_Bool SCIPsdpiSolverIsPrimalInfeasible(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverReadSDP(SCIP_SDPISOLVER *sdpisolver, const char *fname)
struct SCIP_SDPiSolver SCIP_SDPISOLVER
SCIP_RETCODE SCIPsdpiSolverComputePenaltyparam(SCIP_SDPISOLVER *sdpisolver, SCIP_Real maxcoeff, SCIP_Real *penaltyparam)
SCIP_Bool SCIPsdpiSolverIsIterlimExc(SCIP_SDPISOLVER *sdpisolver)
int SCIPsdpiSolverGetInternalStatus(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverLoadAndSolveWithPenalty(SCIP_SDPISOLVER *sdpisolver, SCIP_Real penaltyparam, SCIP_Bool withobj, SCIP_Bool rbound, int nvars, SCIP_Real *obj, SCIP_Real *lb, SCIP_Real *ub, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nremovedblocks, int nlpcons, int noldlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int *rownactivevars, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *start, SCIP_SDPSOLVERSETTING startsettings, SCIP_Real timelimit, SCIP_Bool *feasorig, SCIP_Bool *penaltybound)
enum SCIP_SDPParam SCIP_SDPPARAM
SCIP_Bool SCIPsdpiSolverIsDualUnbounded(SCIP_SDPISOLVER *sdpisolver)
SCIP_Bool SCIPsdpiSolverIsTimelimExc(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverLoadAndSolve(SCIP_SDPISOLVER *sdpisolver, int nvars, SCIP_Real *obj, SCIP_Real *lb, SCIP_Real *ub, int nsdpblocks, int *sdpblocksizes, int *sdpnblockvars, int sdpconstnnonz, int *sdpconstnblocknonz, int **sdpconstrow, int **sdpconstcol, SCIP_Real **sdpconstval, int sdpnnonz, int **sdpnblockvarnonz, int **sdpvar, int ***sdprow, int ***sdpcol, SCIP_Real ***sdpval, int **indchanges, int *nremovedinds, int *blockindchanges, int nremovedblocks, int nlpcons, int noldlpcons, SCIP_Real *lplhs, SCIP_Real *lprhs, int *lprownactivevars, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *start, SCIP_SDPSOLVERSETTING startsettings, SCIP_Real timelimit)
SCIP_Bool SCIPsdpiSolverIsDualFeasible(SCIP_SDPISOLVER *sdpisolver)
SCIP_RETCODE SCIPsdpiSolverGetIterations(SCIP_SDPISOLVER *sdpisolver, int *iterations)
SCIP_Bool SCIPsdpiSolverIsConverged(SCIP_SDPISOLVER *sdpisolver)