SCIP-SDP  2.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sdpisolver_sdpa.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of SCIPSDP - a solving framework for mixed-integer */
4 /* semidefinite programms based on SCIP. */
5 /* */
6 /* Copyright (C) 2011-2013 Discrete Optimization, TU Darmstadt */
7 /* EDOM, FAU Erlangen-Nürnberg */
8 /* 2014-2016 Discrete Optimization, TU Darmstadt */
9 /* */
10 /* */
11 /* This program is free software; you can redistribute it and/or */
12 /* modify it under the terms of the GNU Lesser General Public License */
13 /* as published by the Free Software Foundation; either version 3 */
14 /* of the License, or (at your option) any later version. */
15 /* */
16 /* This program is distributed in the hope that it will be useful, */
17 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
18 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
19 /* GNU Lesser General Public License for more details. */
20 /* */
21 /* You should have received a copy of the GNU Lesser General Public License */
22 /* along with this program; if not, write to the Free Software */
23 /* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
24 /* */
25 /* */
26 /* Based on SCIP - Solving Constraint Integer Programs */
27 /* Copyright (C) 2002-2016 Zuse Institute Berlin */
28 /* SCIP is distributed under the terms of the SCIP Academic Licence, */
29 /* see file COPYING in the SCIP distribution. */
30 /* */
31 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
32 
33 /* #define SCIP_DEBUG*/
34 /* #define SCIP_MORE_DEBUG *//* shows all added nonzero entries */
35 /* #define SCIP_DEBUG_PRINTTOFILE *//* prints each problem inserted into SDPA to the file sdpa.dat-s and the starting point to sdpa.ini-s */
36 
37 /* #define SDPA_RESETPARAMS */ /* this can be used together with an update to the SDPA source code to prevent memory leaks when using SCIP-SDP with SDPA */
38 
45 /* turn off lint warnings for whole file: */
46 /*lint --e{750,788,818}*/
47 
48 #include <assert.h>
49 
50 #include "sdpi/sdpisolver.h"
51 
52 /* turn off warnings for sdpa (doesn't seem to work) */
53 #pragma GCC diagnostic ignored "-Wstrict-prototypes"
54 #include "sdpa_call.h" /* SDPA callable library interface */
55 #pragma GCC diagnostic warning "-Wstrict-prototypes"
56 
57 #include "blockmemshell/memory.h" /* for memory allocation */
58 #include "scip/def.h" /* for SCIP_Real, _Bool, ... */
59 #include "scip/pub_misc.h" /* for sorting */
60 
61 /* local defines */
62 #define EPSILONCHANGE 1
63 #define FEASTOLCHANGE 1
64 #define PENALTYBOUNDTOL 1E-3
67 #define MIN_LAMBDASTAR 1e0
68 #define MAX_LAMBDASTAR 1e8
69 #define LAMBDASTAR_FACTOR 1e0
70 #define LAMBDASTAR_TWOPOINTS TRUE
71 #define LAMBDASTAR_THRESHOLD 1e1
72 #define LAMBDASTAR_LOW 1.5
73 #define LAMBDASTAR_HIGH 1e5
75 #define MIN_PENALTYPARAM 1e5
76 #define MAX_PENALTYPARAM 1e12
77 #define PENALTYPARAM_FACTOR 1e1
78 #define MAX_MAXPENALTYPARAM 1e15
79 #define MAXPENALTYPARAM_FACTOR 1e6
82 #define BMS_CALL(x) do \
83  { \
84  if( NULL == (x) ) \
85  { \
86  SCIPerrorMessage("No memory in function call.\n"); \
87  return SCIP_NOMEMORY; \
88  } \
89  } \
90  while( FALSE )
91 
93 #define CHECK_IF_SOLVED(sdpisolver) do \
94  { \
95  if (!(sdpisolver->solved)) \
96  { \
97  SCIPerrorMessage("Tried to access solution information for SDP %d ahead of solving!\n", sdpisolver->sdpcounter); \
98  return SCIP_LPERROR; \
99  } \
100  } \
101  while( FALSE )
102 
104 #define CHECK_IF_SOLVED_BOOL(sdpisolver) do \
105  { \
106  if (!(sdpisolver->solved)) \
107  { \
108  SCIPerrorMessage("Tried to access solution information for SDP %d ahead of solving!\n", sdpisolver->sdpcounter); \
109  return FALSE; \
110  } \
111  } \
112  while( FALSE )
113 
114 
116 struct SCIP_SDPiSolver
117 {
118  SCIP_MESSAGEHDLR* messagehdlr;
119  BMS_BLKMEM* blkmem;
120  SDPA* sdpa;
121  int nvars;
122  int nactivevars;
123  int* inputtosdpamapper;
126  int* sdpatoinputmapper;
127  SCIP_Real* fixedvarsval;
128  SCIP_Real fixedvarsobjcontr;
129  int nvarbounds;
130  int* varboundpos;
132  SCIP_Bool solved;
133  SCIP_Bool timelimit;
134  int sdpcounter;
135  int niterations;
136  int nsdpcalls;
137  SCIP_Real epsilon;
138  SCIP_Real feastol;
139  SCIP_Real objlimit;
140  SCIP_Bool sdpinfo;
141  SCIP_Bool penalty;
142  SCIP_Bool rbound;
143  SCIP_SDPSOLVERSETTING usedsetting;
144  SCIP_Real lambdastar;
145 };
146 
147 
148 /*
149  * Local Functions
150  */
151 
152 #ifndef NDEBUG
153 
154 static
155 SCIP_Bool isFixed(
156  SCIP_SDPISOLVER* sdpisolver,
157  SCIP_Real lb,
158  SCIP_Real ub
159  )
160 {
161  assert( sdpisolver != NULL );
162  assert( lb < ub + sdpisolver->feastol );
163 
164  return (REALABS(ub-lb) <= sdpisolver->feastol);
165 }
166 #else
167 #define isFixed(sdpisolver,lb,ub) (REALABS(ub-lb) <= sdpisolver->feastol)
168 #endif
169 
170 
171 /*
172  * Miscellaneous Methods
173  */
174 
180 const char* SCIPsdpiSolverGetSolverName(
181  void
182  )
183 {/*lint !e1784*/
184  return "SDPA"; /* version number can only be printed to file/stdout */
185 }
186 
188 const char* SCIPsdpiSolverGetSolverDesc(
189  void
190  )
191 {/*lint !e1784*/
192  return "Primal-dual Interior Point Solver for SDPs developed by K. Fujisawa et al. (sdpa.sourceforge.net)";
193 }
194 
201 void* SCIPsdpiSolverGetSolverPointer(
202  SCIP_SDPISOLVER* sdpisolver
203  )
204 {/*lint !e1784*/
205  assert( sdpisolver != NULL );
206  return (void*) sdpisolver->sdpa;
207 }
208 
212 /*
213  * SDPI Creation and Destruction Methods
214  */
215 
220 SCIP_RETCODE SCIPsdpiSolverCreate(
221  SCIP_SDPISOLVER** sdpisolver,
222  SCIP_MESSAGEHDLR* messagehdlr,
223  BMS_BLKMEM* blkmem
224  )
225 {/*lint !e1784*/
226  assert( sdpisolver != NULL );
227  assert( blkmem != NULL );
228 
229  SCIPdebugMessage("Calling SCIPsdpiCreate \n");
230 
231  BMS_CALL( BMSallocBlockMemory(blkmem, sdpisolver) );
232 
233  (*sdpisolver)->messagehdlr = messagehdlr;
234  (*sdpisolver)->blkmem = blkmem;
235 
236  /* this will be properly initialized then calling solve */
237  (*sdpisolver)->sdpa = NULL;
238 
239  (*sdpisolver)->nvars = 0;
240  (*sdpisolver)->nactivevars = 0;
241  (*sdpisolver)->inputtosdpamapper = NULL;
242  (*sdpisolver)->sdpatoinputmapper = NULL;
243  (*sdpisolver)->fixedvarsval = NULL;
244  (*sdpisolver)->fixedvarsobjcontr = 0.0;
245  (*sdpisolver)->nvarbounds = 0;
246  (*sdpisolver)->varboundpos = NULL;
247  (*sdpisolver)->solved = FALSE;
248  (*sdpisolver)->timelimit = FALSE;
249  (*sdpisolver)->sdpcounter = 0;
250  (*sdpisolver)->niterations = 0;
251  (*sdpisolver)->nsdpcalls = 0;
252 
253  (*sdpisolver)->epsilon = 1e-4;
254  (*sdpisolver)->feastol = 1e-6;
255  (*sdpisolver)->objlimit = SCIPsdpiSolverInfinity(*sdpisolver);
256  (*sdpisolver)->sdpinfo = FALSE;
257  (*sdpisolver)->usedsetting = SCIP_SDPSOLVERSETTING_UNSOLVED;
258 
259  return SCIP_OKAY;
260 }
261 
263 SCIP_RETCODE SCIPsdpiSolverFree(
264  SCIP_SDPISOLVER** sdpisolver
265  )
266 {/*lint !e1784*/
267  assert( sdpisolver != NULL );
268  assert( *sdpisolver != NULL );
269 
270  SCIPdebugMessage("Freeing SDPISolver\n");
271 
272  /* free SDPA object using destructor and free memory via blockmemshell */
273  if ( (*sdpisolver)->sdpa != NULL)
274  delete (*sdpisolver)->sdpa;
275 
276  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->varboundpos, 2 * (*sdpisolver)->nvars); /*lint !e647*/
277 
278  if ( (*sdpisolver)->nvars > 0 )
279  BMSfreeBlockMemoryArray((*sdpisolver)->blkmem, &(*sdpisolver)->inputtosdpamapper, (*sdpisolver)->nvars);/*lint !e737*/
280 
281  if ( (*sdpisolver)->nactivevars > 0 )
282  BMSfreeBlockMemoryArray((*sdpisolver)->blkmem, &(*sdpisolver)->sdpatoinputmapper, (*sdpisolver)->nactivevars);/*lint !e737*/
283 
284  if ( (*sdpisolver)->nvars > (*sdpisolver)->nactivevars )
285  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->fixedvarsval, (*sdpisolver)->nvars - (*sdpisolver)->nactivevars); /*lint !e776*/
286 
287  BMSfreeBlockMemory((*sdpisolver)->blkmem, sdpisolver);
288 
289  return SCIP_OKAY;
290 }
291 
293 SCIP_RETCODE SCIPsdpiSolverIncreaseCounter(
294  SCIP_SDPISOLVER* sdpisolver
295  )
296 {/*lint !e1784*/
297  assert( sdpisolver != NULL );
298 
299  sdpisolver->sdpcounter++;
300 
301  return SCIP_OKAY;
302 }
303 
305 SCIP_RETCODE SCIPsdpiSolverResetCounter(
306  SCIP_SDPISOLVER* sdpisolver
307  )
308 {/*lint !e1784*/
309  assert( sdpisolver != NULL );
310 
311  SCIPdebugMessage("Resetting counter of SDP-Interface from %d to 0.\n", sdpisolver->sdpcounter);
312  sdpisolver->sdpcounter = 0;
313 
314  return SCIP_OKAY;
315 }
316 
320 /*
321  * Solving Methods
322  */
323 
339 SCIP_RETCODE SCIPsdpiSolverLoadAndSolve(
340  SCIP_SDPISOLVER* sdpisolver,
341  int nvars,
342  SCIP_Real* obj,
343  SCIP_Real* lb,
344  SCIP_Real* ub,
345  int nsdpblocks,
346  int* sdpblocksizes,
347  int* sdpnblockvars,
348  int sdpconstnnonz,
349  int* sdpconstnblocknonz,
351  int** sdpconstrow,
352  int** sdpconstcol,
353  SCIP_Real** sdpconstval,
354  int sdpnnonz,
355  int** sdpnblockvarnonz,
357  int** sdpvar,
359  int*** sdprow,
360  int*** sdpcol,
361  SCIP_Real*** sdpval,
362  int** indchanges,
364  int* nremovedinds,
365  int* blockindchanges,
366  int nremovedblocks,
367  int nlpcons,
368  int noldlpcons,
369  SCIP_Real* lplhs,
370  SCIP_Real* lprhs,
371  int* lprownactivevars,
372  int lpnnonz,
373  int* lprow,
374  int* lpcol,
375  SCIP_Real* lpval,
376  SCIP_Real* start,
377  SCIP_SDPSOLVERSETTING startsettings,
379  SCIP_Real timelimit
380  )
381 {/*lint !e1784*/
382  return SCIPsdpiSolverLoadAndSolveWithPenalty(sdpisolver, 0.0, TRUE, FALSE, nvars, obj, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
383  sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
384  nremovedinds, blockindchanges, nremovedblocks, nlpcons, noldlpcons, lplhs, lprhs, lprownactivevars, lpnnonz, lprow, lpcol, lpval, start,
385  startsettings, timelimit, NULL, NULL);
386 }
387 
407 SCIP_RETCODE SCIPsdpiSolverLoadAndSolveWithPenalty(
408  SCIP_SDPISOLVER* sdpisolver,
409  SCIP_Real penaltyparam,
410  SCIP_Bool withobj,
411  SCIP_Bool rbound,
412  int nvars,
413  SCIP_Real* obj,
414  SCIP_Real* lb,
415  SCIP_Real* ub,
416  int nsdpblocks,
417  int* sdpblocksizes,
418  int* sdpnblockvars,
419  int sdpconstnnonz,
420  int* sdpconstnblocknonz,
422  int** sdpconstrow,
423  int** sdpconstcol,
424  SCIP_Real** sdpconstval,
425  int sdpnnonz,
426  int** sdpnblockvarnonz,
428  int** sdpvar,
430  int*** sdprow,
431  int*** sdpcol,
432  SCIP_Real*** sdpval,
433  int** indchanges,
435  int* nremovedinds,
436  int* blockindchanges,
437  int nremovedblocks,
438  int nlpcons,
439  int noldlpcons,
440  SCIP_Real* lplhs,
441  SCIP_Real* lprhs,
442  int* rownactivevars,
443  int lpnnonz,
444  int* lprow,
445  int* lpcol,
446  SCIP_Real* lpval,
447  SCIP_Real* start,
448  SCIP_SDPSOLVERSETTING startsettings,
450  SCIP_Real timelimit,
451  SCIP_Bool* feasorig,
453  SCIP_Bool* penaltybound
455  )
456 {/*lint !e1784*/
457  SCIP_Real* sdpavarbounds;
458  int i;
459  int k;
460  int block;
461  int nfixedvars;
462  bool checkinput; /* should the input be checked for consistency in SDPA ? */
463  int lpconsind;
464  int lastrow;
465  int* rowmapper; /* maps the lhs- and rhs-inequalities of the old LP-cons to their constraint numbers in DSDP */
466  int nlpineqs;
467  int pos;
468  int newpos;
469  int oldnvars;
470 #ifdef SCIP_MORE_DEBUG
471  int ind;
472 #endif
473 
474 #ifdef SCIP_DEBUG
475  char phase_string[15];
476 #endif
477 
478  assert( sdpisolver != NULL );
479  assert( penaltyparam > -1 * sdpisolver->epsilon );
480  assert( penaltyparam < sdpisolver->epsilon || ( feasorig != NULL ) );
481  assert( nvars > 0 );
482  assert( obj != NULL );
483  assert( lb != NULL );
484  assert( ub != NULL );
485  assert( nsdpblocks >= 0 );
486  assert( nsdpblocks == 0 || sdpblocksizes != NULL );
487  assert( nsdpblocks == 0 || sdpnblockvars != NULL );
488  assert( sdpconstnnonz >= 0 );
489  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstnblocknonz != NULL );
490  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstrow != NULL );
491  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstcol != NULL );
492  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstval != NULL );
493  assert( sdpnnonz >= 0 );
494  assert( nsdpblocks == 0 || sdpnblockvarnonz != NULL );
495  assert( nsdpblocks == 0 || sdpvar != NULL );
496  assert( nsdpblocks == 0 || sdprow != NULL );
497  assert( nsdpblocks == 0 || sdpcol != NULL );
498  assert( nsdpblocks == 0 || sdpval != NULL );
499  assert( nsdpblocks == 0 || indchanges != NULL );
500  assert( nsdpblocks == 0 || nremovedinds != NULL );
501  assert( nsdpblocks == 0 || blockindchanges != NULL );
502  assert( 0 <= nremovedblocks && nremovedblocks <= nsdpblocks );
503  assert( nlpcons >= 0 );
504  assert( noldlpcons >= nlpcons );
505  assert( nlpcons == 0 || lplhs != NULL );
506  assert( nlpcons == 0 || lprhs != NULL );
507  assert( nlpcons == 0 || rownactivevars != NULL );
508  assert( lpnnonz >= 0 );
509  assert( nlpcons == 0 || lprow != NULL );
510  assert( nlpcons == 0 || lpcol != NULL );
511  assert( nlpcons == 0 || lpval != NULL );
512 
513  sdpisolver->niterations = 0;
514  sdpisolver->nsdpcalls = 0;
515 
516  /* immediately exit if the time limit is negative */
517  if ( timelimit <= 0.0 )
518  {
519  sdpisolver->solved = FALSE;
520  sdpisolver->timelimit = TRUE;
521  return SCIP_OKAY;
522  }
523  else
524  sdpisolver->timelimit = FALSE;
525 
526 #ifndef SCIP_NDEBUG
527  checkinput = false;
528 #else
529  checkinput = true;
530 #endif
531 
532  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_UNSOLVED;
533 
534  /* only increase the counter if we don't use the penalty formulation to stay in line with the numbers in the general interface (where this is still the
535  * same SDP) */
536  if ( penaltyparam < sdpisolver->epsilon )
537  SCIPdebugMessage("Inserting Data into SDPA for SDP (%d) \n", ++sdpisolver->sdpcounter);
538  else
539  SCIPdebugMessage("Inserting Data again into SDPA for SDP (%d) \n", sdpisolver->sdpcounter);
540 
541  /* set the penalty and rbound flags accordingly */
542  sdpisolver->penalty = (penaltyparam < sdpisolver->epsilon) ? FALSE : TRUE;
543  sdpisolver->rbound = rbound;
544 
545  /* allocate memory for inputtosdpamapper, sdpatoinputmapper and the fixed variable information, for the latter this will
546  * later be shrinked if the needed size is known */
547  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->inputtosdpamapper), sdpisolver->nvars, nvars) );
548  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->sdpatoinputmapper), sdpisolver->nactivevars, nvars) );
549  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->fixedvarsval), sdpisolver->nvars - sdpisolver->nactivevars, nvars) ); /*lint !e776*/
550 
551  oldnvars = sdpisolver->nvars; /* we need to save this to realloc the varboundpos-array if needed */
552  sdpisolver->nvars = nvars;
553  sdpisolver->nactivevars = 0;
554  nfixedvars = 0;
555 
556  /* find the fixed variables */
557  sdpisolver->fixedvarsobjcontr = 0.0;
558  for (i = 0; i < nvars; i++)
559  {
560  if ( isFixed(sdpisolver, lb[i], ub[i]) )
561  {
562  sdpisolver->fixedvarsobjcontr += obj[i] * lb[i]; /* this is the value this fixed variable contributes to the objective */
563  sdpisolver->fixedvarsval[nfixedvars] = lb[i]; /* if lb=ub, than this is the value the variable will have in every solution */
564  nfixedvars++;
565  sdpisolver->inputtosdpamapper[i] = -nfixedvars;
566  SCIPdebugMessage("Fixing variable %d locally to %f for SDP %d in SDPA\n", i, lb[i], sdpisolver->sdpcounter);
567  }
568  else
569  {
570  sdpisolver->sdpatoinputmapper[sdpisolver->nactivevars] = i;
571  sdpisolver->nactivevars++;
572  sdpisolver->inputtosdpamapper[i] = sdpisolver->nactivevars; /* sdpa starts counting at 1, so we do this after increasing nactivevars */
573 #ifdef SCIP_MORE_DEBUG
574  SCIPdebugMessage("Variable %d becomes variable %d for SDP %d in SDPA\n", i, sdpisolver->inputtosdpamapper[i], sdpisolver->sdpcounter);
575 #endif
576  }
577  }
578  assert( sdpisolver->nactivevars + nfixedvars == sdpisolver->nvars );
579 
580  /* if we want to solve without objective, we reset fixedvarsobjcontr */
581  if ( ! withobj )
582  sdpisolver->fixedvarsobjcontr = 0.0;
583 
584  /* shrink the fixedvars and sdpatoinputmapper arrays to the right size */
585  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->fixedvarsval), nvars, nfixedvars) );
586  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->sdpatoinputmapper), nvars, sdpisolver->nactivevars) );
587 
588  /* insert data */
589  if ( sdpisolver->sdpa != 0 )
590  {
591  /* if the SDPA solver has already been created, clear the current problem instance */
592  delete sdpisolver->sdpa;
593  sdpisolver->sdpa = new SDPA();
594  }
595  else
596  sdpisolver->sdpa = new SDPA();
597  assert( sdpisolver->sdpa != 0 );
598 
599  /* initialize settings (this needs to be done before inserting the problem as the initial point depends on the settings) */
600  if ( penaltyparam >= sdpisolver->epsilon || startsettings == SCIP_SDPSOLVERSETTING_STABLE || startsettings == SCIP_SDPSOLVERSETTING_PENALTY )
601  {
602  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_STABLE_BUT_SLOW); /* if we already had problems with this problem, there is no reason to try fast */
603  /* as we want to solve with stable settings, we also update epsilon and the feasibility tolerance, as we skip the default settings, we multpy twice */
604  sdpisolver->sdpa->setParameterEpsilonStar(EPSILONCHANGE * EPSILONCHANGE * sdpisolver->epsilon);
605  sdpisolver->sdpa->setParameterEpsilonDash(FEASTOLCHANGE * FEASTOLCHANGE * sdpisolver->feastol);
606  SCIPdebugMessage("Start solving process with stable settings\n");
607  }
608  else if ( startsettings == SCIP_SDPSOLVERSETTING_UNSOLVED || startsettings == SCIP_SDPSOLVERSETTING_FAST)
609  {
610  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_UNSTABLE_BUT_FAST);
611  sdpisolver->sdpa->setParameterEpsilonStar(sdpisolver->epsilon);
612  sdpisolver->sdpa->setParameterEpsilonDash(sdpisolver->feastol);
613  SCIPdebugMessage("Start solving process with fast settings\n");
614  }
615  else if ( startsettings == SCIP_SDPSOLVERSETTING_MEDIUM )
616  {
617  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_DEFAULT);
618  /* as we want to solve with stable settings, we also update epsilon and the feasibility tolerance, as we skip the default settings, we multpy once */
619  sdpisolver->sdpa->setParameterEpsilonStar(EPSILONCHANGE * sdpisolver->epsilon);
620  sdpisolver->sdpa->setParameterEpsilonDash(FEASTOLCHANGE * sdpisolver->feastol);
621  SCIPdebugMessage("Start solving process with medium settings\n");
622  }
623  else
624  {
625  SCIPdebugMessage("Unknown setting for start-settings: %d!\n", startsettings); \
626  return SCIP_LPERROR;
627  }
628  sdpisolver->sdpa->setParameterLowerBound(-1e20);
629 
630 
631  /* set the objective limit */
632  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
633  sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
634  else
635  sdpisolver->sdpa->setParameterUpperBound(1e8);
636 #ifdef SCIP_MORE_DEBUG
637  sdpisolver->sdpa->printParameters(stdout);
638 #endif
639 
640  /* increase Lambda Star, this seems to help the numerics */
641  sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
642 
643  /* compute number of variable bounds and save them in sdpavarbounds */
644  sdpisolver->nvarbounds = 0;
645  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &sdpavarbounds, 2 * sdpisolver->nactivevars) ); /*lint !e647*//*lint !e530*/
646 
647  if ( sdpisolver->nvars != oldnvars )
648  {
649  if ( sdpisolver->varboundpos == NULL )
650  {
651  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->varboundpos), 2 * sdpisolver->nvars) ); /*lint !e647*/
652  }
653  else
654  {
655  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &(sdpisolver->varboundpos), oldnvars, 2 * sdpisolver->nvars) ); /*lint !e647*/
656  }
657  }
658  assert( sdpisolver->varboundpos != NULL );
659 
660  for (i = 0; i < sdpisolver->nactivevars; i++)
661  {
662  assert( 0 <= sdpisolver->sdpatoinputmapper[i] && sdpisolver->sdpatoinputmapper[i] < nvars );
663  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, lb[sdpisolver->sdpatoinputmapper[i]]))
664  {
665  sdpavarbounds[sdpisolver->nvarbounds] = lb[sdpisolver->sdpatoinputmapper[i]];
666  sdpisolver->varboundpos[sdpisolver->nvarbounds] = -(i + 1); /* negative sign means lower bound, i + 1 because sdpa starts counting from one */
667  (sdpisolver->nvarbounds)++;
668  }
669  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, ub[sdpisolver->sdpatoinputmapper[i]]) )
670  {
671  sdpavarbounds[sdpisolver->nvarbounds] = ub[sdpisolver->sdpatoinputmapper[i]];
672  sdpisolver->varboundpos[sdpisolver->nvarbounds] = +(i + 1); /* positive sign means upper bound, i + 1 because sdpa starts counting from one */
673  (sdpisolver->nvarbounds)++;
674  }
675  }
676 
677  if ( nlpcons > 0 )
678  {
679  /* allocate memory to save which lpconstraints are mapped to which index, entry 2i corresponds to the left hand side of row i, 2i+1 to the rhs */
680  BMS_CALL( BMSallocBlockMemoryArray(sdpisolver->blkmem, &rowmapper, 2*noldlpcons) ); /*lint !e647*//*lint !e530*/
681 
682  /* compute the number of LP constraints after splitting the ranged rows and compute the rowmapper */
683  pos = 1; /* SDPA starts counting the LP-inequalities at one */
684  newpos = 0; /* the position in the lhs and rhs arrays */
685  for (i = 0; i < noldlpcons; i++)
686  {
687  if ( rownactivevars[i] >= 2 )
688  {
689  if ( lplhs[newpos] > - SCIPsdpiSolverInfinity(sdpisolver) )
690  {
691  rowmapper[2*i] = pos; /*lint !e679*/
692  pos++;
693  }
694  else
695  rowmapper[2*i] = -1; /*lint !e679*/
696  if ( lprhs[newpos] < SCIPsdpiSolverInfinity(sdpisolver) )
697  {
698  rowmapper[2*i + 1] = pos; /*lint !e679*/
699  pos++;
700  }
701  else
702  rowmapper[2*i + 1] = -1; /*lint !e679*/
703 
704  newpos++;
705  }
706  else
707  {
708  rowmapper[2*i] = -1; /*lint !e679*/
709  rowmapper[2*i + 1] = -1; /*lint !e679*/
710  }
711  }
712  nlpineqs = pos - 1; /* minus one because we started at one as SDPA wants them numbered one to nlpineqs */
713  assert( nlpineqs <= 2*nlpcons ); /* *2 comes from left- and right-hand-sides */
714  }
715  else
716  nlpineqs = 0;
717 
718  /* if we use a penalty formulation, we need the constraint r >= 0 */
719  if ( penaltyparam >= sdpisolver->epsilon && rbound )
720  sdpisolver->nvarbounds++;
721 
722  if ( sdpisolver->sdpinfo )
723  sdpisolver->sdpa->setDisplay(stdout);
724  else
725  sdpisolver->sdpa->setDisplay(0);
726 
727 #ifdef SCIP_MORE_DEBUG
728  FILE* fpOut = fopen("output.tmp", "w");
729  if ( ! fpOut )
730  exit(-1);
731  sdpisolver->sdpa->setResultFile(fpOut);
732 #endif
733 
734  /* initialize blockstruct */
735  if ( penaltyparam < sdpisolver->epsilon ) /* we initialize this with an exact 0.0 in Solve without penalty */
736  sdpisolver->sdpa->inputConstraintNumber((long long) sdpisolver->nactivevars);/*lint !e747*/
737  else
738  sdpisolver->sdpa->inputConstraintNumber((long long) sdpisolver->nactivevars + 1);/*lint !e747*/ /* the additional variable is r for the penalty formulation */
739 
740  /* if there are any lp-cons/variable-bounds, we get an extra block for those, lastrow - nshifts is the number of lp constraints added */
741  sdpisolver->sdpa->inputBlockNumber((long long) ((nlpineqs + sdpisolver->nvarbounds > 0) ?
742  nsdpblocks - nremovedblocks + 1 : nsdpblocks - nremovedblocks)); /*lint !e834 !e747 !e776*/
743 
744  /* block+1 because SDPA starts counting at 1 */
745  for (block = 0; block < nsdpblocks; block++)
746  {
747  if ( blockindchanges[block] >= 0 )
748  {
749  SCIPdebugMessage("adding block %d to SDPA as block %d with size %d\n",
750  block, block - blockindchanges[block] + 1, sdpblocksizes[block] - nremovedinds[block]); /*lint !e834*/
751  sdpisolver->sdpa->inputBlockSize((long long) block - blockindchanges[block] + 1,/*lint !e747, !e834*/
752  (long long) sdpblocksizes[block] - nremovedinds[block]); /*lint !e834, !e776, !e747*/
753  sdpisolver->sdpa->inputBlockType((long long) block - blockindchanges[block] + 1, SDPA::SDP); /*lint !e834, !e776, !e747*/
754  }
755  }
756  if ( nlpineqs + sdpisolver->nvarbounds > 0 )
757  {
758  /* the last block is the lp block, the size has a negative sign */
759  sdpisolver->sdpa->inputBlockSize((long long) nsdpblocks - nremovedblocks + 1, (long long) -(nlpineqs + sdpisolver->nvarbounds)); /*lint !e834, !e776, !e747*/
760  sdpisolver->sdpa->inputBlockType((long long) nsdpblocks - nremovedblocks + 1, SDPA::LP); /*lint !e834, !e776, !e747*/
761  SCIPdebugMessage("adding LP block to SDPA as block %d with size %d\n", nsdpblocks - nremovedblocks + 1,/*lint !e834*/
762  -(nlpineqs + sdpisolver->nvarbounds)); /*lint !e834*/
763  }
764 
765  sdpisolver->sdpa->initializeUpperTriangleSpace();
766 
767  /* set objective values */
768  for (i = 0; i < sdpisolver->nactivevars; i++)
769  {
770  if ( withobj )
771  {
772  /* insert objective value, SDPA counts from 1 to n instead of 0 to n-1 */
773  sdpisolver->sdpa->inputCVec((long long) i + 1, obj[sdpisolver->sdpatoinputmapper[i]]);/*lint !e747*/
774 #ifdef SCIP_MORE_DEBUG
775  SCIPdebugMessage("inserting objective %f for variable %d which became variable %d in SDPA\n", obj[sdpisolver->sdpatoinputmapper[i]],
776  sdpisolver->sdpatoinputmapper[i], i+1);
777 #endif
778  }
779  if ( penaltyparam >= sdpisolver->epsilon )
780  sdpisolver->sdpa->inputCVec((long long) sdpisolver->nactivevars + 1, penaltyparam);/*lint !e747*/ /* set the objective of the additional var to penaltyparam */
781  }
782 
783  /* if we want to use a starting point we have to tell SDPA to allocate memory for it */
784  if ( start != NULL )
785  sdpisolver->sdpa->setInitPoint(true);
786  else
787  sdpisolver->sdpa->setInitPoint(false);
788 
789  /* start inserting the non-constant SDP-Constraint-Matrices */
790  if ( sdpnnonz > 0 )
791  {
792  int v;
793  int blockvar;
794 
795  assert( nsdpblocks > 0 );
796  assert( sdpnblockvarnonz != NULL );
797  assert( sdpnblockvars != NULL );
798  assert( sdpcol != NULL );
799  assert( sdprow != NULL );
800  assert( sdpval != NULL );
801  assert( sdpvar != NULL );
802  assert( indchanges != NULL );
803  assert( nremovedinds != NULL );
804 
805  for (block = 0; block < nsdpblocks; block++)
806  {
807  /* if the block has no entries, we skip it */
808  if ( blockindchanges[block] == -1 )
809  continue;
810 #ifdef SCIP_MORE_DEBUG
811  SCIPdebugMessage(" -> building block %d, which becomes block %d in SDPA (%d)\n", block, block - blockindchanges[block] + 1,sdpisolver->sdpcounter);
812 #endif
813  /* iterate over all variables in this block */
814  for (blockvar = 0; blockvar < sdpnblockvars[block]; blockvar++)
815  {
816  v = sdpisolver->inputtosdpamapper[sdpvar[block][blockvar]];
817 
818 #ifdef SCIP_MORE_DEBUG
819  SCIPdebugMessage(" -> adding coefficient matrix for variable %d which becomes variable %d in SDPA (%d)\n",
820  sdpvar[block][blockvar], v, sdpisolver->sdpcounter);
821 #endif
822 
823  /* check if the variable is active */
824  if ( v > -1 )
825  {
826  for (k = 0; k < sdpnblockvarnonz[block][blockvar]; k++)
827  {
828  /* rows and cols with active nonzeros should not be removed */
829  assert( indchanges[block][sdprow[block][blockvar][k]] > -1 && indchanges[block][sdpcol[block][blockvar][k]] > -1 );
830 
831  assert( indchanges[block][sdprow[block][blockvar][k]] <= sdprow[block][blockvar][k]);
832  assert( indchanges[block][sdpcol[block][blockvar][k]] <= sdpcol[block][blockvar][k]);
833 
834  assert( 0 <= sdprow[block][blockvar][k] && sdprow[block][blockvar][k] < sdpblocksizes[block] );
835  assert( 0 <= sdpcol[block][blockvar][k] && sdpcol[block][blockvar][k] < sdpblocksizes[block] );
836 
837  /* rows and columns start with one in SDPA, so we have to add 1 to the indices */
838 #ifdef SCIP_MORE_DEBUG
839  SCIPdebugMessage(" -> adding nonzero %g at (%d,%d) (%d)\n",
840  sdpval[block][blockvar][k],
841  sdpcol[block][blockvar][k] - indchanges[block][sdpcol[block][blockvar][k]] + 1,
842  sdprow[block][blockvar][k] - indchanges[block][sdprow[block][blockvar][k]] + 1,
843  sdpisolver->sdpcounter);
844 #endif
845 
846  sdpisolver->sdpa->inputElement((long long) v, (long long) block - blockindchanges[block] + 1, /*lint !e834, !e747*/
847  (long long) sdpcol[block][blockvar][k] - indchanges[block][sdpcol[block][blockvar][k]] + 1, /*lint !e834, !e747*/
848  (long long) sdprow[block][blockvar][k] - indchanges[block][sdprow[block][blockvar][k]] + 1, /*lint !e834, !e747*/
849  sdpval[block][blockvar][k], checkinput); /*lint !e776 !e834*/
850  }
851  }
852  }
853  /* insert the identity matrix if we are using a penalty formulation */
854  if ( penaltyparam >= sdpisolver->epsilon )
855  {
856 #ifdef SCIP_MORE_DEBUG
857  SCIPdebugMessage(" -> adding coefficient matrix for penalty variable r in SDPA (%d)\n", sdpisolver->sdpcounter);
858 #endif
859  for (i = 0; i < sdpblocksizes[block] - nremovedinds[block]; i++)
860  {
861 #ifdef SCIP_MORE_DEBUG
862  SCIPdebugMessage(" -> adding nonzero 1.0 at (%d,%d) (%d)\n", i + 1, i + 1, sdpisolver->sdpcounter);
863 #endif
864 
865  sdpisolver->sdpa->inputElement((long long) sdpisolver->nactivevars + 1, (long long) block - blockindchanges[block] + 1, /*lint !e747, !e776, !e834*/
866  (long long) i + 1, (long long) i + 1, 1.0, checkinput);/*lint !e747*/
867  }
868  }
869  }
870  }
871 
872  /* start inserting the constant matrix */
873  if ( sdpconstnnonz > 0 )
874  {
875  assert( nsdpblocks > 0 );
876  assert( sdpconstnblocknonz!= NULL );
877  assert( sdpconstcol != NULL );
878  assert( sdpconstrow != NULL );
879  assert( sdpconstval != NULL );
880 
881  for (block = 0; block < nsdpblocks; block++)
882  {
883  if ( blockindchanges[block] == -1 )
884  continue;
885 #ifdef SCIP_MORE_DEBUG
886  SCIPdebugMessage(" -> building block %d (%d)\n", block + 1, sdpisolver->sdpcounter);
887 #endif
888  for (k = 0; k < sdpconstnblocknonz[block]; k++)
889  {
890  /* rows and cols with active nonzeros should not be removed */
891  assert( indchanges[block][sdpconstrow[block][k]] > -1 && indchanges[block][sdpconstcol[block][k]] > -1 );
892 
893  assert( indchanges[block][sdpconstrow[block][k]] <= sdpconstrow[block][k]);
894  assert( indchanges[block][sdpconstcol[block][k]] <= sdpconstcol[block][k]);
895 
896  assert (0 <= sdpconstrow[block][k] && sdpconstrow[block][k] < sdpblocksizes[block]);
897  assert (0 <= sdpconstcol[block][k] && sdpconstcol[block][k] < sdpblocksizes[block]);
898 
899  /* rows and columns start with one in SDPA, so we have to add 1 to the indices, the constant matrix is given as variable 0 */
900 #ifdef SCIP_MORE_DEBUG
901  SCIPdebugMessage(" -> adding constant nonzero %g at (%d,%d) (%d)\n", sdpconstval[block][k],
902  sdpconstcol[block][k] - indchanges[block][sdpconstcol[block][k]] + 1,
903  sdpconstrow[block][k] - indchanges[block][sdpconstrow[block][k]] + 1,
904  sdpisolver->sdpcounter);
905 #endif
906  sdpisolver->sdpa->inputElement((long long) 0, (long long) block - blockindchanges[block] + 1, /*lint !e747, !e776, !e834*/
907  (long long) sdpconstcol[block][k] - indchanges[block][sdpconstcol[block][k]] + 1, /*lint !e747, !e776, !e834*/
908  (long long) sdpconstrow[block][k] - indchanges[block][sdpconstrow[block][k]] + 1, /*lint !e747, !e776, !e834*/
909  sdpconstval[block][k], checkinput);
910  }
911  }
912  }
913 
914 #ifdef SCIP_MORE_DEBUG
915  SCIPdebugMessage(" -> building LP-block %d (%d)\n", nsdpblocks - nremovedblocks + 1, sdpisolver->sdpcounter);
916 #endif
917  /* inserting LP nonzeros */
918  lastrow = -1;
919  for (i = 0; i < lpnnonz; i++)
920  {
921  assert( 0 <= lprow[i] && lprow[i] < noldlpcons );
922  assert( 0 <= lpcol[i] && lpcol[i] < nvars );
923  assert( REALABS(lpval[i]) > sdpisolver->feastol );
924 
925  /* if the variable is active and the constraint is more than a bound, we add it */
926  if ( sdpisolver->inputtosdpamapper[lpcol[i]] > 0 )
927  {
928  /* as this is an active variable, there should be at least one in the constraint */
929  assert( rownactivevars[lprow[i]] > 0 );
930  if ( rownactivevars[lprow[i]] > 1 )
931  {
932  if ( lprow[i] > lastrow ) /* we update the lpcons-counter */
933  {
934  lastrow = lprow[i];
935  /* if we use a penalty formulation, add the r * Identity entry */
936  if ( penaltyparam >= sdpisolver->epsilon )
937  {
938  /* check for the lhs-inequality */
939  if ( rowmapper[2*lastrow] > -1 ) /*lint !e679*/
940  {
941 #ifdef SCIP_MORE_DEBUG
942  SCIPdebugMessage(" -> adding nonzero 1.0 at (%d,%d) for penalty variable r in SDPA (%d)\n",
943  rowmapper[2*lastrow], rowmapper[2*lastrow], sdpisolver->sdpcounter);
944 #endif
945  /* LP nonzeros are added as diagonal entries of the last block (coming after the last SDP-block, with
946  * blocks starting at 1, as are rows), the r-variable is variable nactivevars + 1 */
947  sdpisolver->sdpa->inputElement((long long) sdpisolver->nactivevars + 1, (long long) nsdpblocks - nremovedblocks + 1, /*lint !e747, !e776, !e834*/
948  (long long) rowmapper[2*lastrow], (long long) rowmapper[2*lastrow], 1.0, checkinput); /*lint !e679, !e747, !e834*/
949  }
950 
951  /* check for the rhs-inequality */
952  if ( rowmapper[2*lastrow + 1] > -1 ) /*lint !e679*/
953  {
954 #ifdef SCIP_MORE_DEBUG
955  SCIPdebugMessage(" -> adding nonzero 1.0 at (%d,%d) for penalty variable r in SDPA (%d)\n",
956  rowmapper[2*lastrow + 1], rowmapper[2*lastrow + 1], sdpisolver->sdpcounter);
957 #endif
958  /* LP nonzeros are added as diagonal entries of the last block (coming after the last SDP-block, with
959  * blocks starting at 1, as are rows), the r-variable is variable nactivevars + 1 */
960  sdpisolver->sdpa->inputElement((long long) sdpisolver->nactivevars + 1, (long long) nsdpblocks - nremovedblocks + 1, /*lint !e747, !e776, !e834*/
961  (long long) rowmapper[2*lastrow + 1], (long long) rowmapper[2*lastrow + 1], 1.0, checkinput); /*lint !e679, !e747, !e834*/
962  }
963  }
964  }
965  /* add the lp-nonzero to the lhs-inequality if it exists: */
966  if ( rowmapper[2*lastrow] > -1 ) /*lint !e679*/
967  {
968 #ifdef SCIP_MORE_DEBUG
969  SCIPdebugMessage(" -> adding nonzero %g at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
970  lpval[i], rowmapper[2*lastrow], rowmapper[2*lastrow], lpcol[i], sdpisolver->inputtosdpamapper[lpcol[i]], sdpisolver->sdpcounter);
971 #endif
972  /* LP nonzeros are added as diagonal entries of the last block (coming after the last SDP-block, with blocks starting at 1, as are rows) */
973  sdpisolver->sdpa->inputElement((long long) sdpisolver->inputtosdpamapper[lpcol[i]], (long long) nsdpblocks - nremovedblocks + 1, /*lint !e747, !e776, !e834*/
974  (long long) rowmapper[2*lastrow], (long long) rowmapper[2*lastrow], lpval[i], checkinput); /*lint !e679, !e747, !e834*/
975  }
976  /* add the lp-nonzero to the rhs-inequality if it exists: */
977  if ( rowmapper[2*lastrow + 1] > -1 ) /*lint !e679*/
978  {
979 #ifdef SCIP_MORE_DEBUG
980  SCIPdebugMessage(" -> adding nonzero %g at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
981  -1 * lpval[i], rowmapper[2*lastrow + 1], rowmapper[2*lastrow + 1], lpcol[i], sdpisolver->inputtosdpamapper[lpcol[i]], sdpisolver->sdpcounter);
982 #endif
983  /* LP nonzeros are added as diagonal entries of the last block (coming after the last SDP-block, with blocks starting at 1, as are rows),
984  * the -1 comes from the fact that this is a <=-constraint, while SDPA works with >= */
985  sdpisolver->sdpa->inputElement((long long) sdpisolver->inputtosdpamapper[lpcol[i]], (long long) nsdpblocks - nremovedblocks + 1, /*lint !e747, !e776, !e834*/
986  (long long) rowmapper[2*lastrow + 1], (long long) rowmapper[2*lastrow + 1], -1 * lpval[i], checkinput); /*lint !e679, !e747, !e834*/
987  }
988  }
989  }
990  }
991 
992  /* inserting LP left- and right-hand-sides for active constraints */
993  lpconsind = 1; /* this is the same order we used when computing the rowmapper, so we insert at the right positions */
994  for (i = 0; i < nlpcons; i++)
995  {
996  /* check for left-hand side */
997  if ( lplhs[i] > - SCIPsdpiSolverInfinity(sdpisolver) )
998  {
999 #ifdef SCIP_MORE_DEBUG
1000  SCIPdebugMessage(" -> adding lhs %g at (%d,%d) (%d)\n", lplhs[i], lpconsind, lpconsind, sdpisolver->sdpcounter);
1001 #endif
1002  /* LP constraints are added as diagonal entries of the last block, left-hand-side is added as variable zero */
1003  sdpisolver->sdpa->inputElement((long long) 0, (long long) nsdpblocks - nremovedblocks + 1, (long long) lpconsind, /*lint !e747, !e776, !e834*/
1004  (long long) lpconsind, lplhs[i], checkinput);/*lint !e747*/
1005  lpconsind++;
1006  }
1007 
1008  /* check for right-hand side */
1009  if ( lprhs[i] < SCIPsdpiSolverInfinity(sdpisolver) )
1010  {
1011 #ifdef SCIP_MORE_DEBUG
1012  SCIPdebugMessage(" -> adding lhs (originally rhs) %g at (%d,%d) (%d)\n", -1 * lprhs[i], lpconsind, lpconsind, sdpisolver->sdpcounter);
1013 #endif
1014  /* LP constraints are added as diagonal entries of the last block, right-hand side is added as variable zero, the -1 comes from
1015  * the fact, that SDPA uses >=, while the rhs is a <=-constraint */
1016  sdpisolver->sdpa->inputElement((long long) 0, (long long) nsdpblocks - nremovedblocks + 1, (long long) lpconsind, /*lint !e747, !e776, !e834*/
1017  (long long) lpconsind, -1 * lprhs[i], checkinput);/*lint !e747*/
1018  lpconsind++;
1019  }
1020  }
1021 
1022  assert( lpconsind == nlpineqs + 1 ); /* plus one because we started at one as SDPA wants them numbered one to nlpineqs */
1023 
1024  /* print each LP-constraint as one formatted constraint in addition to the single entries inserted into SDPA */
1025 #ifdef SCIP_MORE_DEBUG
1026  lastrow = -1;
1027  ind = -1; /* this is increased once before the first usage */
1028  for (i = 0; i < lpnnonz; i++)
1029  {
1030  /* if the variable is active and the constraint is more than a bound, we added it */
1031  if ( sdpisolver->inputtosdpamapper[lpcol[i]] > 0 )
1032  {
1033  if ( rownactivevars[lprow[i]] > 1 )
1034  {
1035  if ( lprow[i] > lastrow ) /* we finished the old row */
1036  {
1037  if ( lastrow >= 0 )
1038  {
1039  if ( lprhs[ind] < SCIPsdpiSolverInfinity(sdpisolver) )
1040  SCIPmessagePrintInfo(sdpi->messagehdlr, " <= %f\n", lprhs[ind]);
1041  else
1042  SCIPmessagePrintInfo(sdpi->messagehdlr, "\n");
1043  }
1044  lastrow = lprow[i];
1045  ind++;
1046  if ( lplhs[ind] > - SCIPsdpiSolverInfinity(sdpisolver) )
1047  SCIPmessagePrintInfo(sdpi->messagehdlr, "%f <= ", lplhs[ind]);
1048  }
1049  SCIPmessagePrintInfo(sdpi->messagehdlr, "+ %f <x%d> ", lpval[i], lpcol[i]);
1050  }
1051  }
1052  }
1053  if ( lastrow >= 0 )
1054  {
1055  if ( lprhs[ind] < SCIPsdpiSolverInfinity(sdpisolver) )
1056  SCIPmessagePrintInfo(sdpi->messagehdlr, " <= %f\n", lprhs[ind]);
1057  else
1058  SCIPmessagePrintInfo(sdpi->messagehdlr, "\n");
1059  }
1060  assert( ind == nlpcons - 1 ); /* -1 because we start indexing at zero and do not increase after the last row */
1061 #endif
1062 
1063  /* free the memory for the rowmapper */
1064  if ( nlpcons > 0 )
1065  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &rowmapper, 2 * noldlpcons);/*lint !e647, !e737*/
1066 
1067  /* insert variable bounds, these are also added as LP-constraints and therefore diagonal entries of the LP block
1068  * if we work with the penalty formulation, we get an extra entry for r >= 0, but this we will add afterwards */
1069  for (i = 0; i < ((penaltyparam < sdpisolver->epsilon) || (! rbound) ? sdpisolver->nvarbounds : sdpisolver->nvarbounds - 1); i++)
1070  {
1071  assert( 0 < abs(sdpisolver->varboundpos[i]) && abs(sdpisolver->varboundpos[i] <= sdpisolver->nactivevars) ); /* the indices are already those for SDPA */
1072 
1073  /* for lower bound */
1074  if ( sdpisolver->varboundpos[i] < 0 )
1075  {
1076  /* add it as an lp-constraint for this variable (- because we saved -n for the lower bound), at the position
1077  * (nactivelpcons + 1) + varbound-index, because we have >= the variable has coefficient +1 */
1078  sdpisolver->sdpa->inputElement((long long) -sdpisolver->varboundpos[i], (long long) nsdpblocks - nremovedblocks + 1, /*lint !e747, !e776, !e834*/
1079  (long long) nlpineqs + 1 + i, (long long) nlpineqs + 1 + i, 1.0, checkinput);/*lint !e747*/
1080 
1081  if ( REALABS(sdpavarbounds[i]) > sdpisolver->epsilon )
1082  {
1083  /* the bound is added as the rhs and therefore variable zero */
1084 #ifdef SCIP_MORE_DEBUG
1085  SCIPdebugMessage(" -> adding lower bound %g at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1086  sdpavarbounds[i], nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpatoinputmapper[-sdpisolver->varboundpos[i] - 1],
1087  -sdpisolver->varboundpos[i], sdpisolver->sdpcounter);
1088 #endif
1089  sdpisolver->sdpa->inputElement((long long) 0, (long long) nsdpblocks - nremovedblocks + 1, (long long) nlpineqs + 1 + i, /*lint !e747, !e776, !e834*/
1090  (long long) nlpineqs + 1 + i, sdpavarbounds[i], checkinput);/*lint !e747*/
1091  }
1092  else
1093  {
1094  /* as the bound is zero, we don't need to add a right hand side */
1095 #ifdef SCIP_MORE_DEBUG
1096  SCIPdebugMessage(" -> adding lower bound 0 at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1097  nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpatoinputmapper[-sdpisolver->varboundpos[i] - 1],
1098  -sdpisolver->varboundpos[i], sdpisolver->sdpcounter);
1099 #endif
1100  }
1101  }
1102  else
1103  {
1104  /* this is an upper bound */
1105 
1106  /* add it as an lp-constraint for this variable, at the position nactivelpcons + varbound-index, because we have >= but we
1107  * want <= for the upper bound, we have to multiply by -1 and therefore the variable has coefficient -1 */
1108  sdpisolver->sdpa->inputElement((long long) sdpisolver->varboundpos[i], (long long) nsdpblocks - nremovedblocks + 1, /*lint !e747, !e776, !e834*/
1109  (long long) nlpineqs + 1 + i, (long long) nlpineqs + 1 + i, -1.0, checkinput);/*lint !e747*/
1110 
1111  if ( REALABS(sdpavarbounds[i]) > sdpisolver->epsilon )
1112  {
1113  /* the bound is added as the rhs and therefore variable zero, we multiply by -1 for <= */
1114 #ifdef SCIP_MORE_DEBUG
1115  SCIPdebugMessage(" -> adding upper bound %g at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1116  sdpavarbounds[i], nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpatoinputmapper[sdpisolver->varboundpos[i] - 1],
1117  sdpisolver->varboundpos[i], sdpisolver->sdpcounter);
1118 #endif
1119  sdpisolver->sdpa->inputElement((long long) 0, (long long) nsdpblocks - nremovedblocks + 1, (long long) nlpineqs + 1 + i, /*lint !e747, !e776, !e834*/
1120  (long long) nlpineqs + 1 + i, -sdpavarbounds[i], checkinput);/*lint !e747*/
1121  }
1122  else
1123  {
1124  /* as the bound is zero, we don't need to add a right hand side */
1125 #ifdef SCIP_MORE_DEBUG
1126  SCIPdebugMessage(" -> adding upper bound 0 at (%d,%d) for variable %d which became variable %d in SDPA (%d)\n",
1127  0, nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpatoinputmapper[sdpisolver->varboundpos[i] - 1],
1128  sdpisolver->varboundpos[i]);
1129 #endif
1130  }
1131  }
1132  }
1133 
1134  if ( penaltyparam >= sdpisolver->epsilon && rbound )
1135  {
1136  /* we add the variable bound r >= 0 */
1137  sdpisolver->sdpa->inputElement((long long) sdpisolver->nactivevars + 1, (long long) nsdpblocks - nremovedblocks + 1, /*lint !e747, !e776, !e834*/
1138  (long long) nlpineqs + 1 + i, (long long) nlpineqs + 1 + i, 1.0, checkinput);/*lint !e747*/
1139 #ifdef SCIP_MORE_DEBUG
1140  SCIPdebugMessage(" -> adding lower bound r >= 0 at (%d,%d) in SDPA (%d)\n", nlpineqs + 1 + i, nlpineqs + 1 + i, sdpisolver->sdpcounter);
1141 #endif
1142  }
1143 
1144  /* free the arrays used for counting and saving variable bounds and LP-right-hand-sides */
1145  BMSfreeBlockMemoryArray(sdpisolver->blkmem, &sdpavarbounds, 2 * sdpisolver->nactivevars); /*lint !e647, !e737*/
1146 
1147  /* transform the matrices to a more efficient form */
1148  sdpisolver->sdpa->initializeUpperTriangle();
1149  sdpisolver->sdpa->initializeSolve();
1150 
1151  /* set the starting solution */
1152  if ( start != NULL && penaltyparam < sdpisolver->epsilon )
1153  {
1154  SCIPdebugMessage("Starting with a previous solution is not yet tested for the interface, only x-vector is given, not y and Z");
1155  for (i = 1; i <= sdpisolver->nactivevars; i++) /* we iterate over the variables in sdpa */
1156  sdpisolver->sdpa->inputInitXVec((long long) i, start[sdpisolver->sdpatoinputmapper[i] - 1]);/*lint !e747*/
1157  }
1158  else if ( penaltyparam >= sdpisolver->epsilon )
1159  SCIPdebugMessage("Skipping insertion of starting point, as this is not yet supported for penalty formulation.\n");
1160 
1161 #ifdef SCIP_DEBUG_PRINTTOFILE
1162  /* if necessary, dump input data and initial point */
1163  sdpisolver->sdpa->writeInputSparse(const_cast<char*>("sdpa.dat-s"), const_cast<char*>("%+8.3e"));
1164  sdpisolver->sdpa->writeInitSparse(const_cast<char*>("sdpa.ini-s"), const_cast<char*>("%+8.3e"));
1165 #endif
1166 
1167  SCIPdebugMessage("Calling SDPA solve (SDP: %d)\n", sdpisolver->sdpcounter);
1168  sdpisolver->sdpa->solve();
1169  sdpisolver->solved = TRUE;
1170 
1171  /* update number of SDP-iterations and -calls */
1172  sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1173  sdpisolver->nsdpcalls += 1;
1174 
1175 #ifdef SCIP_DEBUG
1176  /* print the phase value , i.e. whether solving was successfull */
1177  sdpisolver->sdpa->getPhaseString((char*)phase_string);
1178  SCIPdebugMessage("SDPA solving finished with status %s (primal and dual here are switched in contrast to our formulation)\n", phase_string);
1179 #endif
1180 
1181  /* remember settings */
1182  if ( SCIPsdpiSolverIsAcceptable(sdpisolver) && penaltyparam < sdpisolver->epsilon )
1183  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_FAST;
1184  else if ( penaltyparam >= sdpisolver->epsilon )
1185  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_PENALTY;
1186 
1187  /* check whether problem has been stably solved, if it wasn't and we didn't yet use the default parametersettings (for the penalty formulation we do so), try
1188  * again with more stable parameters */
1189  if ( ! SCIPsdpiSolverIsAcceptable(sdpisolver) && penaltyparam < sdpisolver->epsilon &&
1190  (startsettings == SCIP_SDPSOLVERSETTING_UNSOLVED || startsettings == SCIP_SDPSOLVERSETTING_FAST) )
1191  {
1192  SCIPdebugMessage("Numerical troubles -- solving SDP %d again ...\n", sdpisolver->sdpcounter);
1193 
1194  /* initialize settings */
1195  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_DEFAULT);
1196  sdpisolver->sdpa->setParameterEpsilonStar(EPSILONCHANGE * sdpisolver->epsilon);
1197  sdpisolver->sdpa->setParameterEpsilonDash(FEASTOLCHANGE * sdpisolver->feastol);
1198  sdpisolver->sdpa->setParameterLowerBound(-1e20);
1199  /* set the objective limit */
1200  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
1201  sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
1202  else
1203  sdpisolver->sdpa->setParameterUpperBound(1e8);
1204 
1205  /* increase Lambda Star, this seems to help the numerics */
1206  sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
1207 
1208 #ifdef SCIP_MORE_DEBUG
1209  sdpisolver->sdpa->printParameters(stdout);
1210 #endif
1211  sdpisolver->sdpa->setInitPoint(false);
1212 #ifdef SDPA_RESETPARAMS
1213  sdpisolver->sdpa->resetParameters();
1214 #else
1215  sdpisolver->sdpa->initializeSolve();
1216 #endif
1217  sdpisolver->sdpa->solve();
1218  sdpisolver->solved = TRUE;
1219 
1220  /* update number of SDP-iterations and -calls */
1221  sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1222  sdpisolver->nsdpcalls += 1;
1223 
1224  /* remember setting */
1225  if ( SCIPsdpiSolverIsAcceptable(sdpisolver) )
1226  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_MEDIUM;
1227 
1228 #ifdef SCIP_DEBUG
1229  /* print the phase value , i.e. whether solving was successfull */
1230  sdpisolver->sdpa->getPhaseString((char*)phase_string);
1231  SCIPdebugMessage("SDPA solving finished with status %s (primal and dual here are switched in contrast to our formulation)\n", phase_string);
1232 #endif
1233 
1234  /* if we still didn't converge, and did not yet use the stable settings, set the parameters even more conservativly */
1235  if ( (! SCIPsdpiSolverIsAcceptable(sdpisolver)) && penaltyparam < sdpisolver->epsilon &&
1236  (startsettings == SCIP_SDPSOLVERSETTING_UNSOLVED || startsettings == SCIP_SDPSOLVERSETTING_FAST || startsettings == SCIP_SDPSOLVERSETTING_MEDIUM) )
1237  {
1238  SCIPdebugMessage("Numerical troubles -- solving SDP %d again^2 ...\n", sdpisolver->sdpcounter);
1239 
1240  /* initialize settings */
1241  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_STABLE_BUT_SLOW);
1242  sdpisolver->sdpa->setParameterEpsilonStar(EPSILONCHANGE * EPSILONCHANGE * sdpisolver->epsilon);
1243  sdpisolver->sdpa->setParameterEpsilonDash(FEASTOLCHANGE * FEASTOLCHANGE * sdpisolver->feastol);
1244  sdpisolver->sdpa->setParameterLowerBound(-1e20);
1245  /* set the objective limit */
1246  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
1247  sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
1248  else
1249  sdpisolver->sdpa->setParameterUpperBound(1e8);
1250 
1251  /* increase Lambda Star, this seems to help the numerics */
1252  sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
1253 
1254 #ifdef SCIP_MORE_DEBUG
1255  sdpisolver->sdpa->printParameters(stdout);
1256 #endif
1257  sdpisolver->sdpa->setInitPoint(false);
1258 #ifdef SDPA_RESETPARAMS
1259  sdpisolver->sdpa->resetParameters();
1260 #else
1261  sdpisolver->sdpa->initializeSolve();
1262 #endif
1263  sdpisolver->sdpa->solve();
1264  sdpisolver->solved = TRUE;
1265 
1266  /* update number of SDP-iterations and -calls */
1267  sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1268  sdpisolver->nsdpcalls += 1;
1269 
1270  /* remember setting */
1271  if ( SCIPsdpiSolverIsAcceptable(sdpisolver) )
1272  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_STABLE;
1273 
1274 #ifdef SCIP_DEBUG
1275  /* print the phase value , i.e. whether solving was successfull */
1276  sdpisolver->sdpa->getPhaseString((char*)phase_string);
1277  SCIPdebugMessage("SDPA solving finished with status %s (primal and dual here are switched in constrast to our formulation)\n", phase_string);
1278 #endif
1279  }
1280  }
1281 
1282 #ifdef SCIP_MORE_DEBUG
1283  (void) fclose(fpOut);
1284 #endif
1285 
1286  /* if we solved a penalty formulation, check if the solution is feasible for the original problem (which is the case iff r < feastol) */
1287  if ( penaltyparam >= sdpisolver->epsilon )
1288  {
1289  SCIP_Real* sdpasol;
1290 
1291  /* in the second case we have r as an additional variable */
1292  assert( (sdpisolver->nactivevars + 1 == sdpisolver->sdpa->getConstraintNumber()) ); /*lint !e776*/
1293 
1294  sdpasol = sdpisolver->sdpa->getResultXVec();
1295 
1296  /* we get r as the last variable in SDPA */
1297  *feasorig = (sdpasol[sdpisolver->nactivevars] < sdpisolver->feastol); /*lint !e413*/
1298 
1299  /* if r > 0 or we are in debug mode, also check the primal bound */
1300 #ifdef NDEBUG
1301  if ( ! *feasorig && penaltybound != NULL )
1302  {
1303 #endif
1304  SCIP_Real* X;
1305  int b;
1306  int nblockssdpa;
1307  int nrow;
1308  SCIP_Real trace = 0.0;
1309 
1310  SCIPdebugMessage("Solution not feasible in original problem, r = %f\n", sdpasol[sdpisolver->nactivevars]);
1311 
1312  /* compute Tr(X) */
1313 
1314  /* iterate over all blocks (SDPA starts counting at one and includes the LP block) */
1315  nblockssdpa = (int) sdpisolver->sdpa->getBlockNumber();
1316  for (b = 1; b <= nblockssdpa; b++)
1317  {
1318  /* get the block from SDPA */
1319  X = sdpisolver->sdpa->getResultYMat((long long) b);/*lint !e747*/
1320  nrow = (int) sdpisolver->sdpa->getBlockSize((long long) b);/*lint !e747*/
1321  assert( nrow >= 0 );
1322 
1323  /* if it is the LP-block, we omit the variable bounds as the penalty variable is not added to them */
1324  if ( sdpisolver->sdpa->getBlockType((long long) b) == SDPA::LP )/*lint !e747*/
1325  {
1326  /* iterate over all diagonal entries (until we reach the varbound part), adding them to the trace */
1327  for (i = 0; i < nrow - sdpisolver->nvarbounds; i++)
1328  trace += X[i]; /* get entry (i+1,i+1) for the diagonal matrix X */
1329  }
1330  else
1331  {
1332  /* iterate over all diagonal entries and add them to the trace */
1333  for (i = 0; i < nrow; i++)
1334  trace += X[i + i*nrow];/*lint !e679*/ /* get entry (i+1,i+1) in X */
1335  }
1336  }
1337 
1338  /* if the relative gap is smaller than the tolerance, we return equality */
1339  if ( (penaltyparam - trace) / penaltyparam < PENALTYBOUNDTOL )/*lint !e414*/
1340  {
1341  if ( penaltybound != NULL )
1342  *penaltybound = TRUE;
1343  SCIPdebugMessage("Tr(X) = %f == %f = Gamma, penalty formulation not exact, Gamma should be increased or problem is infeasible\n",
1344  trace, penaltyparam);
1345  }
1346  else if ( penaltybound != NULL )
1347  *penaltybound = FALSE;
1348 
1349 #ifdef NDEBUG
1350  }
1351 #endif
1352  }
1353 
1354  return SCIP_OKAY;
1355 }
1356 
1362 /*
1363  * Solution Information Methods
1364  */
1365 
1370 SCIP_Bool SCIPsdpiSolverWasSolved(
1371  SCIP_SDPISOLVER* sdpisolver
1372  )
1373 {/*lint !e1784*/
1374  assert( sdpisolver != NULL );
1375  return sdpisolver->solved;
1376 }
1377 
1384 SCIP_Bool SCIPsdpiSolverFeasibilityKnown(
1385  SCIP_SDPISOLVER* sdpisolver
1386  )
1387 {/*lint !e1784*/
1388  SDPA::PhaseType phasetype;
1389 
1390  assert( sdpisolver != NULL );
1391  assert( sdpisolver->sdpa != NULL);
1392  CHECK_IF_SOLVED_BOOL( sdpisolver );
1393 
1394  phasetype = sdpisolver->sdpa->getPhaseValue();
1395 
1396  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1397  return FALSE;
1398 
1399  return TRUE;
1400 }
1401 
1403 SCIP_RETCODE SCIPsdpiSolverGetSolFeasibility(
1404  SCIP_SDPISOLVER* sdpisolver,
1405  SCIP_Bool* primalfeasible,
1406  SCIP_Bool* dualfeasible
1407  )
1408 {/*lint !e1784*/
1409  SDPA::PhaseType phasetype;
1410 
1411  assert( sdpisolver != NULL );
1412  assert( sdpisolver->sdpa != NULL );
1413  assert( primalfeasible != NULL );
1414  assert( dualfeasible != NULL );
1415  CHECK_IF_SOLVED( sdpisolver );
1416 
1417  phasetype = sdpisolver->sdpa->getPhaseValue();
1418 
1419  switch ( phasetype )/*lint --e{788}*/
1420  {
1421  case SDPA::pdOPT:
1422  *primalfeasible = TRUE;
1423  *dualfeasible = TRUE;
1424  break;
1425  case SDPA::pdFEAS:
1426  *primalfeasible = TRUE;
1427  *dualfeasible = TRUE;
1428  break;
1429  case SDPA::pFEAS_dINF:
1430  *primalfeasible = TRUE;
1431  *dualfeasible = FALSE;
1432  break;
1433  case SDPA::pINF_dFEAS:
1434  *primalfeasible = FALSE;
1435  *dualfeasible = TRUE;
1436  break;
1437  case SDPA::pUNBD:
1438  *primalfeasible = TRUE;
1439  *dualfeasible = FALSE;
1440  SCIPdebugMessage("SDPA stopped because dual objective became smaller than lower bound\n");
1441  break;
1442  case SDPA::dUNBD:
1443  *primalfeasible = FALSE;
1444  *dualfeasible = TRUE;
1445  SCIPdebugMessage("SDPA stopped because primal objective became bigger than upper bound\n");
1446  break;
1447  default: /* contains noInfo, pFeas, dFeas, pdInf */
1448  SCIPerrorMessage("SDPA doesn't know if primal and dual solutions are feasible\n");
1449  return SCIP_LPERROR;
1450  }
1451 
1452  return SCIP_OKAY;
1453 }
1454 
1457 SCIP_Bool SCIPsdpiSolverIsPrimalUnbounded(
1458  SCIP_SDPISOLVER* sdpisolver
1459  )
1460 {/*lint !e1784*/
1461  SDPA::PhaseType phasetype;
1462 
1463  assert( sdpisolver != NULL );
1464  assert( sdpisolver->sdpa != NULL );
1465  CHECK_IF_SOLVED_BOOL( sdpisolver );
1466 
1467  phasetype = sdpisolver->sdpa->getPhaseValue();
1468 
1469  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::pdINF )
1470  {
1471  SCIPdebugMessage("SDPA doesn't know if primal problem is unbounded");
1472  return FALSE;
1473  }
1474  else if ( phasetype == SDPA::pFEAS_dINF )
1475  return TRUE;
1476  else if ( phasetype == SDPA::pUNBD )
1477  {
1478  SCIPdebugMessage("SDPA was stopped because primal objective became bigger than upper bound");
1479  return TRUE;
1480  }
1481 
1482  return FALSE;
1483 }
1484 
1487 SCIP_Bool SCIPsdpiSolverIsPrimalInfeasible(
1488  SCIP_SDPISOLVER* sdpisolver
1489  )
1490 {/*lint !e1784*/
1491  SDPA::PhaseType phasetype;
1492 
1493  assert( sdpisolver != NULL );
1494  assert( sdpisolver->sdpa != NULL );
1495  CHECK_IF_SOLVED_BOOL( sdpisolver );
1496 
1497  phasetype = sdpisolver->sdpa->getPhaseValue();
1498 
1499  if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1500  {
1501  SCIPdebugMessage("SDPA doesn't know if primal problem is infeasible");
1502  return FALSE;
1503  }
1504  else if ( phasetype == SDPA::pINF_dFEAS )
1505  return TRUE;
1506  else if ( phasetype == SDPA::dUNBD )
1507  {
1508  SCIPdebugMessage("SDPA was stopped because dual objective became smaller than lower bound");
1509  return TRUE;
1510  }
1511 
1512  return FALSE;
1513 }
1514 
1517 SCIP_Bool SCIPsdpiSolverIsPrimalFeasible(
1518  SCIP_SDPISOLVER* sdpisolver
1519  )
1520 {/*lint !e1784*/
1521  SDPA::PhaseType phasetype;
1522 
1523  assert( sdpisolver != NULL );
1524  assert( sdpisolver->sdpa != NULL );
1525  CHECK_IF_SOLVED_BOOL( sdpisolver );
1526 
1527  phasetype = sdpisolver->sdpa->getPhaseValue();
1528 
1529  if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1530  {
1531  SCIPdebugMessage("SDPA doesn't know if primal problem is feasible");
1532  return FALSE;
1533  }
1534  else if ( phasetype == SDPA::pFEAS_dINF || phasetype == SDPA::pdOPT || phasetype == SDPA::pFEAS || phasetype == SDPA::pdFEAS )
1535  return TRUE;
1536  else if ( phasetype == SDPA::dUNBD )
1537  {
1538  SCIPdebugMessage("SDPA was stopped because dual objective became smaller than lower bound");
1539  return TRUE;
1540  }
1541 
1542  return FALSE;
1543 }
1544 
1547 SCIP_Bool SCIPsdpiSolverIsDualUnbounded(
1548  SCIP_SDPISOLVER* sdpisolver
1549  )
1550 {/*lint !e1784*/
1551  SDPA::PhaseType phasetype;
1552 
1553  assert( sdpisolver != NULL );
1554  assert( sdpisolver->sdpa != NULL);
1555  CHECK_IF_SOLVED_BOOL( sdpisolver );
1556 
1557  phasetype = sdpisolver->sdpa->getPhaseValue();
1558 
1559  if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1560  {
1561  SCIPdebugMessage("SDPA doesn't know if dual problem is unbounded");
1562  return FALSE;
1563  }
1564  else if ( phasetype == SDPA::pINF_dFEAS )
1565  return TRUE;
1566  else if ( phasetype == SDPA::dUNBD )
1567  {
1568  SCIPdebugMessage("SDPA was stopped because dual objective became smaller than lower bound");
1569  return TRUE;
1570  }
1571 
1572  return FALSE;
1573 }
1574 
1577 SCIP_Bool SCIPsdpiSolverIsDualInfeasible(
1578  SCIP_SDPISOLVER* sdpisolver
1579  )
1580 {/*lint !e1784*/
1581  SDPA::PhaseType phasetype;
1582 
1583  assert( sdpisolver != NULL );
1584  assert( sdpisolver->sdpa != NULL);
1585  CHECK_IF_SOLVED_BOOL( sdpisolver );
1586 
1587  phasetype = sdpisolver->sdpa->getPhaseValue();
1588 
1589  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::pdINF )
1590  {
1591  SCIPdebugMessage("SDPA doesn't know if dual problem is infeasible");
1592  return FALSE;
1593  }
1594  else if ( phasetype == SDPA::pFEAS_dINF )
1595  return TRUE;
1596  else if ( phasetype == SDPA::pUNBD )
1597  {
1598  SCIPdebugMessage("SDPA was stopped because primal objective became bigger than upper bound");
1599  return TRUE;
1600  }
1601 
1602  return FALSE;
1603 }
1604 
1607 SCIP_Bool SCIPsdpiSolverIsDualFeasible(
1608  SCIP_SDPISOLVER* sdpisolver
1609  )
1610 {/*lint !e1784*/
1611  SDPA::PhaseType phasetype;
1612 
1613  assert( sdpisolver != NULL );
1614  assert( sdpisolver->sdpa != NULL);
1615  CHECK_IF_SOLVED_BOOL( sdpisolver );
1616 
1617  phasetype = sdpisolver->sdpa->getPhaseValue();
1618 
1619  if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1620  {
1621  SCIPdebugMessage("SDPA doesn't know if primal problem is feasible");
1622  return FALSE;
1623  }
1624  else if ( phasetype == SDPA::pINF_dFEAS || phasetype == SDPA::pdOPT || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
1625  return TRUE;
1626  else if ( phasetype == SDPA::dUNBD )
1627  {
1628  SCIPdebugMessage("SDPA was stopped because dual objective became smaller than lower bound");
1629  return TRUE;
1630  }
1631 
1632  return FALSE;
1633 }
1634 
1636 SCIP_Bool SCIPsdpiSolverIsConverged(
1637  SCIP_SDPISOLVER* sdpisolver
1638  )
1639 {/*lint !e1784*/
1640  SDPA::PhaseType phasetype;
1641 
1642  assert( sdpisolver != NULL );
1643  assert( sdpisolver->sdpa != NULL );
1644  CHECK_IF_SOLVED_BOOL( sdpisolver );
1645 
1646  phasetype = sdpisolver->sdpa->getPhaseValue();
1647 
1648  if ( phasetype == SDPA::pdOPT )
1649  return TRUE;
1650 
1651  return FALSE;
1652 }
1653 
1655 SCIP_Bool SCIPsdpiSolverIsObjlimExc(
1656  SCIP_SDPISOLVER* sdpisolver
1657  )
1658 {/*lint !e1784*/
1659  SDPA::PhaseType phasetype;
1660 
1661  assert( sdpisolver != NULL );
1662  assert( sdpisolver->sdpa != NULL );
1663  CHECK_IF_SOLVED_BOOL( sdpisolver );
1664 
1665  phasetype = sdpisolver->sdpa->getPhaseValue();
1666 
1667  if ( phasetype == SDPA::pUNBD )
1668  return TRUE;
1669 
1670  return FALSE;
1671 }
1672 
1674 SCIP_Bool SCIPsdpiSolverIsIterlimExc(
1675  SCIP_SDPISOLVER* sdpisolver
1676  )
1677 {/*lint !e1784*/
1678  SDPA::PhaseType phasetype;
1679 
1680  assert( sdpisolver != NULL );
1681  assert( sdpisolver->sdpa != NULL);
1682  CHECK_IF_SOLVED_BOOL( sdpisolver );
1683 
1684  phasetype = sdpisolver->sdpa->getPhaseValue();
1685 
1686  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
1687  {
1688  if ( sdpisolver->sdpa->getParameterMaxIteration() == sdpisolver->sdpa->getIteration() )
1689  return TRUE;
1690  }
1691 
1692  return FALSE;
1693 }
1694 
1696 SCIP_Bool SCIPsdpiSolverIsTimelimExc(
1697  SCIP_SDPISOLVER* sdpisolver
1698  )
1699 {/*lint !e1784*/
1700  assert( sdpisolver != NULL );
1701 
1702  return sdpisolver->timelimit;
1703 }
1704 
1715 int SCIPsdpiSolverGetInternalStatus(
1716  SCIP_SDPISOLVER* sdpisolver
1717  )
1718 {/*lint !e1784*/
1719  SDPA::PhaseType phasetype;
1720 
1721  assert( sdpisolver != NULL );
1722  assert( sdpisolver->sdpa != NULL );
1723 
1724  if ( sdpisolver->sdpa == NULL || (! sdpisolver->solved) )
1725  return -1;
1726 
1727  phasetype = sdpisolver->sdpa->getPhaseValue();
1728 
1729  if ( phasetype == SDPA::pdOPT || phasetype == SDPA::pFEAS_dINF || phasetype == SDPA::pINF_dFEAS )
1730  return 0;
1731  if ( phasetype == SDPA::pdINF )
1732  return 1;
1733  if ( phasetype == SDPA::pUNBD)
1734  return 3;
1735  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
1736  return 4;
1737  else /* should include dUNBD */
1738  return 7;
1739 }
1740 
1742 SCIP_Bool SCIPsdpiSolverIsOptimal(
1743  SCIP_SDPISOLVER* sdpisolver
1744  )
1745 {/*lint !e1784*/
1746  SDPA::PhaseType phasetype;
1747 
1748  assert( sdpisolver != NULL );
1749  assert( sdpisolver->sdpa != NULL );
1750  CHECK_IF_SOLVED_BOOL( sdpisolver );
1751 
1752  phasetype = sdpisolver->sdpa->getPhaseValue();
1753 
1754  if ( phasetype == SDPA::pdOPT )
1755  return TRUE;
1756 
1757  return FALSE;
1758 }
1759 
1762 SCIP_Bool SCIPsdpiSolverIsAcceptable(
1763  SCIP_SDPISOLVER* sdpisolver
1764  )
1765 {/*lint !e1784*/
1766  SDPA::PhaseType phasetype;
1767 
1768  assert( sdpisolver != NULL );
1769  assert( sdpisolver->sdpa != NULL );
1770 
1771  if ( sdpisolver->timelimit )
1772  return FALSE;
1773 
1774  CHECK_IF_SOLVED_BOOL( sdpisolver );
1775 
1776  phasetype = sdpisolver->sdpa->getPhaseValue();
1777 
1778  /* we are happy if we converged, or we reached the objective limit (pUNBD) or we could show that our problem is
1779  * infeasible [except for numerics], or unbounded */
1780  if ( SCIPsdpiSolverIsConverged(sdpisolver) || phasetype == SDPA::pUNBD || phasetype == SDPA::pINF_dFEAS || phasetype == SDPA::pFEAS_dINF )
1781  return TRUE;
1782 
1783  return FALSE;
1784 }
1785 
1787 SCIP_RETCODE SCIPsdpiSolverIgnoreInstability(
1788  SCIP_SDPISOLVER* sdpisolver,
1789  SCIP_Bool* success
1790  )
1791 {/*lint !e1784*/
1792  SCIPdebugMessage("Not implemented yet\n");
1793 
1794  return SCIP_LPERROR;
1795 }/*lint !e715*/
1796 
1798 SCIP_RETCODE SCIPsdpiSolverGetObjval(
1799  SCIP_SDPISOLVER* sdpisolver,
1800  SCIP_Real* objval
1801  )
1802 {/*lint !e1784*/
1803  assert( sdpisolver != NULL );
1804  assert( sdpisolver->sdpa != NULL );
1805  assert( objval != NULL );
1806  CHECK_IF_SOLVED( sdpisolver );
1807 
1808  *objval = sdpisolver->sdpa->getPrimalObj();
1809 
1810 #ifndef NDEBUG
1811  SCIP_Real primalval = sdpisolver->sdpa->getDualObj();
1812  SCIP_Real gap = (REALABS(*objval - primalval) / (0.5 * (REALABS(primalval) + REALABS(*objval)))); /* duality gap used in SDPA */
1813  if ( gap > sdpisolver->epsilon )
1814  SCIPdebugMessage("Attention: got objective value (before adding values of fixed variables) of %f in SCIPsdpiSolverGetObjval, "
1815  "but primal objective is %f with duality gap %f!\n", *objval, primalval, gap );
1816 #endif
1817 
1818  /* as we didn't add the fixed (lb = ub) variables to sdpa, we have to add their contributions to the objective by hand */
1819  *objval += sdpisolver->fixedvarsobjcontr;
1820 
1821  return SCIP_OKAY;
1822 }
1823 
1828 SCIP_RETCODE SCIPsdpiSolverGetSol(
1829  SCIP_SDPISOLVER* sdpisolver,
1830  SCIP_Real* objval,
1831  SCIP_Real* dualsol,
1832  int* dualsollength
1834  )
1835 {/*lint !e1784*/
1836  SCIP_Real* sdpasol;
1837  int v;
1838 
1839  assert( sdpisolver != NULL );
1840  assert( sdpisolver->sdpa != NULL );
1841  assert( dualsollength != NULL );
1842  CHECK_IF_SOLVED( sdpisolver );
1843 
1844  if ( objval != NULL )
1845  {
1846  *objval = sdpisolver->sdpa->getPrimalObj();
1847 
1848 #ifndef NDEBUG
1849  SCIP_Real primalval = sdpisolver->sdpa->getDualObj();
1850  SCIP_Real gap = (REALABS(*objval - primalval) / (0.5 * (REALABS(primalval) + REALABS(*objval)))); /* duality gap used in SDPA */
1851  if ( gap > sdpisolver->epsilon )
1852  {
1853  SCIPdebugMessage("Attention: got objective value (before adding values of fixed variables) of %f in SCIPsdpiSolverGetSol, "
1854  "but primal objective is %f with duality gap %f!\n", *objval, primalval, gap );
1855  }
1856 #endif
1857 
1858  /* as we didn't add the fixed (lb = ub) variables to sdpa, we have to add their contributions to the objective by hand */
1859  *objval += sdpisolver->fixedvarsobjcontr;
1860  }
1861 
1862  if ( *dualsollength > 0 )
1863  {
1864  assert( dualsol != NULL );
1865  if ( *dualsollength < sdpisolver->nvars )
1866  {
1867  SCIPdebugMessage("The given array in SCIPsdpiSolverGetSol only had length %d, but %d was needed", *dualsollength, sdpisolver->nvars);
1868  *dualsollength = sdpisolver->nvars;
1869 
1870  return SCIP_OKAY;
1871  }
1872 
1873  /* get the solution from sdpa */
1874  assert( (sdpisolver->penalty && sdpisolver->nactivevars + 1 == sdpisolver->sdpa->getConstraintNumber()) || /*lint !e776*/
1875  sdpisolver->nactivevars == sdpisolver->sdpa->getConstraintNumber() ); /* in the second case we have r as an additional variable */
1876  sdpasol = sdpisolver->sdpa->getResultXVec();
1877  /* insert the entries into dualsol, for non-fixed vars we copy those from sdpa, the rest are the saved entries from inserting (they equal lb=ub) */
1878  for (v = 0; v < sdpisolver->nvars; v++)
1879  {
1880  if ( sdpisolver->inputtosdpamapper[v] > 0 )
1881  {
1882  /* minus one because the inputtosdpamapper gives the sdpa indices which start at one, but the array starts at 0 */
1883  dualsol[v] = sdpasol[sdpisolver->inputtosdpamapper[v] - 1];
1884  }
1885  else
1886  {
1887  /* this is the value that was saved when inserting, as this variable has lb=ub */
1888  assert( -sdpisolver->inputtosdpamapper[v] <= sdpisolver->nvars - sdpisolver->nactivevars );
1889  dualsol[v] = sdpisolver->fixedvarsval[(-1 * sdpisolver->inputtosdpamapper[v]) - 1]; /*lint !e679*/
1890  }
1891  }
1892  }
1893  return SCIP_OKAY;
1894 }
1895 
1903 SCIP_RETCODE SCIPsdpiSolverGetPrimalBoundVars(
1904  SCIP_SDPISOLVER* sdpisolver,
1905  SCIP_Real* lbvars,
1906  SCIP_Real* ubvars,
1907  int* arraylength
1909  )
1910 {/*lint !e1784*/
1911  int i;
1912  SCIP_Real* X; /* block of primal solution matrix corresponding to the LP-part */
1913  int lpblockind;
1914  int nlpcons;
1915 
1916  assert( sdpisolver != NULL );
1917  assert( sdpisolver->sdpa != NULL );
1918  assert( lbvars != NULL );
1919  assert( ubvars != NULL );
1920  assert( arraylength != NULL );
1921  assert( *arraylength >= 0 );
1922  CHECK_IF_SOLVED( sdpisolver );
1923 
1924  /* check if the arrays are long enough */
1925  if ( *arraylength < sdpisolver->nvars )
1926  {
1927  *arraylength = sdpisolver->nvars;
1928  SCIPdebugMessage("Insufficient length of array in SCIPsdpiSolverGetPrimalBoundVars (gave %d, needed %d)\n", *arraylength, sdpisolver->nvars);
1929  return SCIP_OKAY;
1930  }
1931 
1932  /* initialize the return-arrays with zero */
1933  for (i = 0; i < sdpisolver->nvars; i++)
1934  {
1935  lbvars[i] = 0.0;
1936  ubvars[i] = 0.0;
1937  }
1938 
1939  /* if no variable bounds were added, we return the zero vector (we do this separately, because in this case there might be no LP-block) */
1940  if ( sdpisolver->nvarbounds == 0 )
1941  {
1942  SCIPdebugMessage("Asked for PrimalBoundVars, but there were no variable bounds in sdpa, returning zero vector !");
1943  return SCIP_OKAY;
1944  }
1945 
1946  /* get the block of primal solution matrix corresponding to the LP-part from sdpa */
1947  lpblockind = (int) sdpisolver->sdpa->getBlockNumber(); /* the LP block is the last one and sdpa counts from one */
1948  assert( sdpisolver->sdpa->getBlockType((long long) lpblockind) == SDPA::LP );/*lint !e747*/
1949  nlpcons = (int) sdpisolver->sdpa->getBlockSize((long long) lpblockind);/*lint !e747*/
1950  assert( nlpcons >= 0 );
1951 
1952  X = sdpisolver->sdpa->getResultYMat((long long) lpblockind);/*lint !e747*/
1953 
1954  /* iterate over all variable bounds and insert the corresponding primal variables in the right positions of the return-arrays */
1955  assert( sdpisolver->nvarbounds <= 2 * sdpisolver->nvars || (sdpisolver->nvarbounds <= 2 * sdpisolver->nvars + 1 && sdpisolver->penalty ) );
1956  /* if we solved a penalty formulation, the last variable bound belongs to the penalty variable, which we aren't interested in here */
1957  for (i = 0; i < ((sdpisolver->penalty) ? sdpisolver->nvarbounds - 1 : sdpisolver->nvarbounds); i++)
1958  {
1959  if ( sdpisolver->varboundpos[i] < 0 )
1960  {
1961  /* this is a lower bound */
1962  /* the last nvarbounds entries correspond to the varbounds */
1963  lbvars[sdpisolver->sdpatoinputmapper[-1 * sdpisolver->varboundpos[i] -1]] = X[nlpcons - sdpisolver->nvarbounds + i]; /*lint !e679, !e834 */
1964  }
1965  else
1966  {
1967  /* this is an upper bound */
1968  /* the last nvarbounds entries correspond to the varbounds */
1969  ubvars[sdpisolver->sdpatoinputmapper[+1 * sdpisolver->varboundpos[i] - 1]] = X[nlpcons - sdpisolver->nvarbounds + i]; /*lint !e679, !e834 */
1970  }
1971  }
1972 
1973  return SCIP_OKAY;
1974 }
1975 
1977 SCIP_RETCODE SCIPsdpiSolverGetIterations(
1978  SCIP_SDPISOLVER* sdpisolver,
1979  int* iterations
1980  )
1981 {/*lint !e1784*/
1982  assert( sdpisolver != NULL );
1983  assert( sdpisolver->sdpa != NULL );
1984  assert( iterations != NULL );
1985 
1986  *iterations = sdpisolver->niterations;
1987 
1988  return SCIP_OKAY;
1989 }
1990 
1992 SCIP_RETCODE SCIPsdpiSolverGetSdpCalls(
1993  SCIP_SDPISOLVER* sdpisolver,
1994  int* calls
1995  )
1996 {/*lint !e1784*/
1997  assert( sdpisolver != NULL );
1998  assert( sdpisolver->sdpa != NULL );
1999  assert( calls != NULL );
2000 
2001  *calls = sdpisolver->nsdpcalls;
2002 
2003  return SCIP_OKAY;
2004 }
2005 
2007 SCIP_RETCODE SCIPsdpiSolverSettingsUsed(
2008  SCIP_SDPISOLVER* sdpisolver,
2009  SCIP_SDPSOLVERSETTING* usedsetting
2010  )
2011 {/*lint !e1784*/
2012  assert( sdpisolver != NULL );
2013  assert( usedsetting != NULL );
2014  CHECK_IF_SOLVED(sdpisolver);
2015 
2016  *usedsetting = sdpisolver->usedsetting;
2017 
2018  return SCIP_OKAY;
2019 }
2020 
2026 /*
2027  * Numerical Methods
2028  */
2029 
2034 SCIP_Real SCIPsdpiSolverInfinity(
2035  SCIP_SDPISOLVER* sdpisolver
2036  )
2037 {/*lint !e1784*/
2038  return 1E+20; /* default infinity from SCIP */
2039 }/*lint !e715*/
2040 
2042 SCIP_Bool SCIPsdpiSolverIsInfinity(
2043  SCIP_SDPISOLVER* sdpisolver,
2044  SCIP_Real val
2045  )
2046 {/*lint !e1784*/
2047  return ( val <= -SCIPsdpiSolverInfinity(sdpisolver) || val >= SCIPsdpiSolverInfinity(sdpisolver) );
2048 }
2049 
2051 SCIP_RETCODE SCIPsdpiSolverGetRealpar(
2052  SCIP_SDPISOLVER* sdpisolver,
2053  SCIP_SDPPARAM type,
2054  SCIP_Real* dval
2055  )
2056 {/*lint !e1784*/
2057  assert( sdpisolver != NULL );
2058  assert( dval != NULL );
2059 
2060  switch( type )/*lint --e{788}*/
2061  {
2062  case SCIP_SDPPAR_EPSILON:
2063  *dval = sdpisolver->epsilon;
2064  break;
2065  case SCIP_SDPPAR_FEASTOL:
2066  *dval = sdpisolver->feastol;
2067  break;
2068  case SCIP_SDPPAR_PENALTYPARAM:
2069  *dval = 0.0;
2070  SCIPdebugMessage("Parameter SCIP_SDPPAR_PENALTYPARAM not used by SDPA"); /* this parameter is only used by DSDP */
2071  break;
2072  case SCIP_SDPPAR_OBJLIMIT:
2073  *dval = sdpisolver->objlimit;
2074  break;
2075  case SCIP_SDPPAR_LAMBDASTAR:
2076  *dval = sdpisolver->lambdastar;
2077  break;
2078  default:
2079  return SCIP_PARAMETERUNKNOWN;
2080  }
2081 
2082  return SCIP_OKAY;
2083 }
2084 
2086 SCIP_RETCODE SCIPsdpiSolverSetRealpar(
2087  SCIP_SDPISOLVER* sdpisolver,
2088  SCIP_SDPPARAM type,
2089  SCIP_Real dval
2090  )
2091 {/*lint !e1784*/
2092  assert( sdpisolver != NULL );
2093 
2094  switch( type )/*lint --e{788}*/
2095  {
2096  case SCIP_SDPPAR_EPSILON:
2097  sdpisolver->epsilon = dval;
2098  SCIPdebugMessage("Setting sdpisolver epsilon to %f.\n", dval);
2099  break;
2100  case SCIP_SDPPAR_FEASTOL:
2101  sdpisolver->feastol = dval;
2102  SCIPdebugMessage("Setting sdpisolver feastol to %f.\n", dval);
2103  break;
2104  case SCIP_SDPPAR_PENALTYPARAM:
2105  SCIPdebugMessage("Parameter SCIP_SDPPAR_PENALTYPARAM not used by SDPA"); /* this parameter is only used by DSDP */
2106  break;
2107  case SCIP_SDPPAR_OBJLIMIT:
2108  SCIPdebugMessage("Setting sdpisolver objlimit to %f.\n", dval);
2109  sdpisolver->objlimit = dval;
2110  break;
2111  case SCIP_SDPPAR_LAMBDASTAR:
2112  SCIPdebugMessage("Setting sdpisolver lambdastar parameter to %f.\n", dval);
2113  sdpisolver->lambdastar = dval;
2114  break;
2115  default:
2116  return SCIP_PARAMETERUNKNOWN;
2117  }
2118 
2119  return SCIP_OKAY;
2120 }
2121 
2123 SCIP_RETCODE SCIPsdpiSolverGetIntpar(
2124  SCIP_SDPISOLVER* sdpisolver,
2125  SCIP_SDPPARAM type,
2126  int* ival
2127  )
2128 {/*lint !e1784*/
2129  assert( sdpisolver != NULL );
2130 
2131  switch( type )/*lint --e{788}*/
2132  {
2133  case SCIP_SDPPAR_SDPINFO:
2134  *ival = (int) sdpisolver->sdpinfo;
2135  SCIPdebugMessage("Getting sdpisolver information output (%d).\n", *ival);
2136  break;
2137  default:
2138  return SCIP_PARAMETERUNKNOWN;
2139  }
2140 
2141  return SCIP_OKAY;
2142 }
2143 
2145 SCIP_RETCODE SCIPsdpiSolverSetIntpar(
2146  SCIP_SDPISOLVER* sdpisolver,
2147  SCIP_SDPPARAM type,
2148  int ival
2149  )
2150 {/*lint !e1784*/
2151  assert( sdpisolver != NULL );
2152 
2153  switch( type )/*lint --e{788}*/
2154  {
2155  case SCIP_SDPPAR_SDPINFO:
2156  sdpisolver->sdpinfo = (SCIP_Bool) ival;
2157  SCIPdebugMessage("Setting sdpisolver information output (%d).\n", ival);
2158  break;
2159  default:
2160  return SCIP_PARAMETERUNKNOWN;
2161  }
2162 
2163  return SCIP_OKAY;
2164 }
2165 
2167 SCIP_RETCODE SCIPsdpiSolverComputeLambdastar(
2168  SCIP_SDPISOLVER* sdpisolver,
2169  SCIP_Real maxguess
2170  )
2171 {/*lint !e1784*/
2172  SCIP_Real compval;
2173 
2174  assert( sdpisolver != NULL );
2175 
2176  /* we set the value to min{max{MIN_LAMBDASTAR, LAMBDASTAR_FACTOR * MAX_GUESS}, MAX_LAMBDASTAR}, where MAX_GUESS is the maximum of the guesses
2177  * of the SDP-Blocks, if the define LAMBDASTAR_TWOPOINTS is set, we instead choose either LAMBDASTAR_LOW or HIGH depending on LAMBASTAR_THRESHOLD */
2178 
2179 #ifdef LAMBDASTAR_TWOPOINTS
2180  if ( maxguess < LAMBDASTAR_THRESHOLD )
2181  compval = LAMBDASTAR_LOW;
2182  else
2183  compval = LAMBDASTAR_HIGH;
2184 #else
2185  compval = LAMBDASTAR_FACTOR * maxguess;
2186 #endif
2187 
2188  if ( compval < MIN_LAMBDASTAR )
2189  {
2190  sdpisolver->lambdastar = MIN_LAMBDASTAR;
2191  SCIPdebugMessage("Setting lambdastar to %f.\n", MIN_LAMBDASTAR);
2192  }
2193  else if ( compval > MAX_LAMBDASTAR )
2194  {
2195  sdpisolver->lambdastar = MAX_LAMBDASTAR;
2196  SCIPdebugMessage("Setting lambdastar to %f.\n", MAX_LAMBDASTAR);
2197  }
2198  else
2199  {
2200  sdpisolver->lambdastar = compval;
2201  SCIPdebugMessage("Setting lambdastar to %f.\n", compval);
2202  }
2203 
2204  return SCIP_OKAY;
2205 }
2206 
2208 SCIP_RETCODE SCIPsdpiSolverComputePenaltyparam(
2209  SCIP_SDPISOLVER* sdpisolver,
2210  SCIP_Real maxcoeff,
2211  SCIP_Real* penaltyparam
2212  )
2213 {/*lint !e1784*/
2214  SCIP_Real compval;
2215 
2216  assert( sdpisolver != NULL );
2217  assert( penaltyparam != NULL );
2218 
2219  compval = PENALTYPARAM_FACTOR * maxcoeff;
2220 
2221  if ( compval < MIN_PENALTYPARAM )
2222  {
2223  SCIPdebugMessage("Setting penaltyparameter to %f.\n", MIN_PENALTYPARAM);
2224  *penaltyparam = MIN_PENALTYPARAM;
2225  }
2226  else if ( compval > MAX_PENALTYPARAM )
2227  {
2228  SCIPdebugMessage("Setting penaltyparameter to %f.\n", MAX_PENALTYPARAM);
2229  *penaltyparam = MAX_PENALTYPARAM;
2230  }
2231  else
2232  {
2233  SCIPdebugMessage("Setting penaltyparameter to %f.\n", compval);
2234  *penaltyparam = compval;
2235  }
2236  return SCIP_OKAY;
2237 }
2238 
2240 SCIP_RETCODE SCIPsdpiSolverComputeMaxPenaltyparam(
2241  SCIP_SDPISOLVER* sdpisolver,
2242  SCIP_Real penaltyparam,
2243  SCIP_Real* maxpenaltyparam
2244  )
2245 {/*lint !e1784*/
2246  SCIP_Real compval;
2247 
2248  assert( sdpisolver != NULL );
2249  assert( maxpenaltyparam != NULL );
2250 
2251  compval = penaltyparam * MAXPENALTYPARAM_FACTOR;
2252 
2253  if ( compval < MAX_MAXPENALTYPARAM )
2254  {
2255  *maxpenaltyparam = compval;
2256  SCIPdebugMessage("Setting maximum penaltyparameter to %f.\n", compval);
2257  }
2258  else
2259  {
2260  *maxpenaltyparam = MAX_MAXPENALTYPARAM;
2261  SCIPdebugMessage("Setting penaltyparameter to %f.\n", MAX_MAXPENALTYPARAM);
2262  }
2263 
2264  return SCIP_OKAY;
2265 }
2266 
2272 /*
2273  * File Interface Methods
2274  */
2275 
2280 SCIP_RETCODE SCIPsdpiSolverReadSDP(
2281  SCIP_SDPISOLVER* sdpisolver,
2282  const char* fname
2283  )
2284 {/*lint !e1784*/
2285  SCIPdebugMessage("Not implemented yet\n");
2286  return SCIP_LPERROR;
2287 }/*lint !e715*/
2288 
2290 SCIP_RETCODE SCIPsdpiSolverWriteSDP(
2291  SCIP_SDPISOLVER* sdpisolver,
2292  const char* fname
2293  )
2294 {/*lint !e1784*/
2295  assert( fname != NULL );
2296 
2297  sdpisolver->sdpa->writeInputSparse(const_cast<char*>(fname), const_cast<char*>("%8.3f"));
2298 
2299  return SCIP_OKAY;
2300 }
2301 
SCIP_SDPSOLVERSETTING_MEDIUM
Definition: type_sdpi.h:71
SCIPsdpiSolverGetSolverName
const char * SCIPsdpiSolverGetSolverName(void)
Definition: sdpisolver_sdpa.cpp:181
MIN_PENALTYPARAM
#define MIN_PENALTYPARAM
Definition: sdpisolver_sdpa.cpp:76
FEASTOLCHANGE
#define FEASTOLCHANGE
Definition: sdpisolver_sdpa.cpp:63
SCIP_SDPSOLVERSETTING_UNSOLVED
Definition: type_sdpi.h:68
SCIPsdpiSolverSetRealpar
SCIP_RETCODE SCIPsdpiSolverSetRealpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, SCIP_Real dval)
Definition: sdpisolver_sdpa.cpp:2087
SCIPsdpiSolverFree
SCIP_RETCODE SCIPsdpiSolverFree(SCIP_SDPISOLVER **sdpisolver)
Definition: sdpisolver_sdpa.cpp:264
SCIPsdpiSolverInfinity
SCIP_Real SCIPsdpiSolverInfinity(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2035
SCIPsdpiSolverGetSolverPointer
void * SCIPsdpiSolverGetSolverPointer(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:202
SCIP_SDPPAR_EPSILON
Definition: type_sdpi.h:53
LAMBDASTAR_LOW
#define LAMBDASTAR_LOW
Definition: sdpisolver_sdpa.cpp:73
LAMBDASTAR_HIGH
#define LAMBDASTAR_HIGH
Definition: sdpisolver_sdpa.cpp:74
SCIPsdpiSolverIgnoreInstability
SCIP_RETCODE SCIPsdpiSolverIgnoreInstability(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *success)
Definition: sdpisolver_sdpa.cpp:1788
SCIP_SDPSOLVERSETTING
enum SCIP_SDPSolverSetting SCIP_SDPSOLVERSETTING
Definition: type_sdpi.h:74
SCIPsdpiSolverIsPrimalUnbounded
SCIP_Bool SCIPsdpiSolverIsPrimalUnbounded(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1458
MAX_LAMBDASTAR
#define MAX_LAMBDASTAR
Definition: sdpisolver_sdpa.cpp:69
SCIPsdpiSolverGetObjval
SCIP_RETCODE SCIPsdpiSolverGetObjval(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *objval)
Definition: sdpisolver_sdpa.cpp:1799
PENALTYPARAM_FACTOR
#define PENALTYPARAM_FACTOR
Definition: sdpisolver_sdpa.cpp:78
SCIPsdpiSolverGetSolverDesc
const char * SCIPsdpiSolverGetSolverDesc(void)
Definition: sdpisolver_sdpa.cpp:189
SCIP_SDPSOLVERSETTING_FAST
Definition: type_sdpi.h:70
SCIPsdpiSolverIsInfinity
SCIP_Bool SCIPsdpiSolverIsInfinity(SCIP_SDPISOLVER *sdpisolver, SCIP_Real val)
Definition: sdpisolver_sdpa.cpp:2043
sdpisolver.h
interface methods for specific SDP-solvers
SCIP_SDPSOLVERSETTING_PENALTY
Definition: type_sdpi.h:69
SCIPsdpiSolverResetCounter
SCIP_RETCODE SCIPsdpiSolverResetCounter(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:306
SCIPsdpiSolverGetRealpar
SCIP_RETCODE SCIPsdpiSolverGetRealpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, SCIP_Real *dval)
Definition: sdpisolver_sdpa.cpp:2052
SCIPsdpiSolverIsObjlimExc
SCIP_Bool SCIPsdpiSolverIsObjlimExc(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1656
EPSILONCHANGE
#define EPSILONCHANGE
Definition: sdpisolver_sdpa.cpp:62
SCIPsdpiSolverGetSol
SCIP_RETCODE SCIPsdpiSolverGetSol(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *objval, SCIP_Real *dualsol, int *dualsollength)
Definition: sdpisolver_sdpa.cpp:1829
MAXPENALTYPARAM_FACTOR
#define MAXPENALTYPARAM_FACTOR
Definition: sdpisolver_sdpa.cpp:80
SCIP_SDPPAR_PENALTYPARAM
Definition: type_sdpi.h:59
SCIPsdpiSolverFeasibilityKnown
SCIP_Bool SCIPsdpiSolverFeasibilityKnown(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1385
MAX_MAXPENALTYPARAM
#define MAX_MAXPENALTYPARAM
Definition: sdpisolver_sdpa.cpp:79
SCIP_SDPPAR_OBJLIMIT
Definition: type_sdpi.h:55
SCIPsdpiSolverComputeLambdastar
SCIP_RETCODE SCIPsdpiSolverComputeLambdastar(SCIP_SDPISOLVER *sdpisolver, SCIP_Real maxguess)
Definition: sdpisolver_sdpa.cpp:2168
SCIPsdpiSolverCreate
SCIP_RETCODE SCIPsdpiSolverCreate(SCIP_SDPISOLVER **sdpisolver, SCIP_MESSAGEHDLR *messagehdlr, BMS_BLKMEM *blkmem)
Definition: sdpisolver_sdpa.cpp:221
SCIPsdpiSolverIsPrimalFeasible
SCIP_Bool SCIPsdpiSolverIsPrimalFeasible(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1518
LAMBDASTAR_THRESHOLD
#define LAMBDASTAR_THRESHOLD
Definition: sdpisolver_sdpa.cpp:72
SCIPsdpiSolverSettingsUsed
SCIP_RETCODE SCIPsdpiSolverSettingsUsed(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPSOLVERSETTING *usedsetting)
Definition: sdpisolver_sdpa.cpp:2008
SCIPsdpiSolverGetIntpar
SCIP_RETCODE SCIPsdpiSolverGetIntpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, int *ival)
Definition: sdpisolver_sdpa.cpp:2124
SCIPsdpiSolverWasSolved
SCIP_Bool SCIPsdpiSolverWasSolved(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1371
LAMBDASTAR_FACTOR
#define LAMBDASTAR_FACTOR
Definition: sdpisolver_sdpa.cpp:70
CHECK_IF_SOLVED_BOOL
#define CHECK_IF_SOLVED_BOOL(sdpisolver)
Definition: sdpisolver_sdpa.cpp:105
isFixed
static SCIP_Bool isFixed(SCIP_SDPISOLVER *sdpisolver, SCIP_Real lb, SCIP_Real ub)
Definition: sdpisolver_sdpa.cpp:156
SCIP_SDPPAR_SDPINFO
Definition: type_sdpi.h:57
MIN_LAMBDASTAR
#define MIN_LAMBDASTAR
Definition: sdpisolver_sdpa.cpp:68
SCIP_SDPPAR_LAMBDASTAR
Definition: type_sdpi.h:61
SCIPsdpiSolverIncreaseCounter
SCIP_RETCODE SCIPsdpiSolverIncreaseCounter(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:294
SCIPsdpiSolverGetSolFeasibility
SCIP_RETCODE SCIPsdpiSolverGetSolFeasibility(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *primalfeasible, SCIP_Bool *dualfeasible)
Definition: sdpisolver_sdpa.cpp:1404
PENALTYBOUNDTOL
#define PENALTYBOUNDTOL
Definition: sdpisolver_sdpa.cpp:64
SCIPsdpiSolverWriteSDP
SCIP_RETCODE SCIPsdpiSolverWriteSDP(SCIP_SDPISOLVER *sdpisolver, const char *fname)
Definition: sdpisolver_sdpa.cpp:2291
SCIP_SDPSOLVERSETTING_STABLE
Definition: type_sdpi.h:72
SCIP_SDPPAR_FEASTOL
Definition: type_sdpi.h:54
SCIPsdpiSolverIsAcceptable
SCIP_Bool SCIPsdpiSolverIsAcceptable(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1763
SCIPsdpiSolverIsOptimal
SCIP_Bool SCIPsdpiSolverIsOptimal(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1743
CHECK_IF_SOLVED
#define CHECK_IF_SOLVED(sdpisolver)
Definition: sdpisolver_sdpa.cpp:94
SCIPsdpiSolverGetSdpCalls
SCIP_RETCODE SCIPsdpiSolverGetSdpCalls(SCIP_SDPISOLVER *sdpisolver, int *calls)
Definition: sdpisolver_sdpa.cpp:1993
SCIPsdpiSolverComputeMaxPenaltyparam
SCIP_RETCODE SCIPsdpiSolverComputeMaxPenaltyparam(SCIP_SDPISOLVER *sdpisolver, SCIP_Real penaltyparam, SCIP_Real *maxpenaltyparam)
Definition: sdpisolver_sdpa.cpp:2241
SCIPsdpiSolverGetPrimalBoundVars
SCIP_RETCODE SCIPsdpiSolverGetPrimalBoundVars(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *lbvars, SCIP_Real *ubvars, int *arraylength)
Definition: sdpisolver_sdpa.cpp:1904
SCIPsdpiSolverIsDualInfeasible
SCIP_Bool SCIPsdpiSolverIsDualInfeasible(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1578
SCIPsdpiSolverSetIntpar
SCIP_RETCODE SCIPsdpiSolverSetIntpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, int ival)
Definition: sdpisolver_sdpa.cpp:2146
SCIPsdpiSolverIsPrimalInfeasible
SCIP_Bool SCIPsdpiSolverIsPrimalInfeasible(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1488
SCIPsdpiSolverReadSDP
SCIP_RETCODE SCIPsdpiSolverReadSDP(SCIP_SDPISOLVER *sdpisolver, const char *fname)
Definition: sdpisolver_sdpa.cpp:2281
SCIP_SDPISOLVER
struct SCIP_SDPiSolver SCIP_SDPISOLVER
Definition: sdpisolver.h:70
SCIPsdpiSolverComputePenaltyparam
SCIP_RETCODE SCIPsdpiSolverComputePenaltyparam(SCIP_SDPISOLVER *sdpisolver, SCIP_Real maxcoeff, SCIP_Real *penaltyparam)
Definition: sdpisolver_sdpa.cpp:2209
SCIPsdpiSolverIsIterlimExc
SCIP_Bool SCIPsdpiSolverIsIterlimExc(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1675
SCIPsdpiSolverGetInternalStatus
int SCIPsdpiSolverGetInternalStatus(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1716
SCIPsdpiSolverLoadAndSolveWithPenalty
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)
Definition: sdpisolver_sdpa.cpp:408
SCIP_SDPPARAM
enum SCIP_SDPParam SCIP_SDPPARAM
Definition: type_sdpi.h:63
SCIPsdpiSolverIsDualUnbounded
SCIP_Bool SCIPsdpiSolverIsDualUnbounded(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1548
SCIPsdpiSolverIsTimelimExc
SCIP_Bool SCIPsdpiSolverIsTimelimExc(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1697
SCIPsdpiSolverLoadAndSolve
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)
Definition: sdpisolver_sdpa.cpp:340
BMS_CALL
#define BMS_CALL(x)
Definition: sdpisolver_sdpa.cpp:83
SCIPsdpiSolverIsDualFeasible
SCIP_Bool SCIPsdpiSolverIsDualFeasible(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1608
SCIPsdpiSolverGetIterations
SCIP_RETCODE SCIPsdpiSolverGetIterations(SCIP_SDPISOLVER *sdpisolver, int *iterations)
Definition: sdpisolver_sdpa.cpp:1978
SCIPsdpiSolverIsConverged
SCIP_Bool SCIPsdpiSolverIsConverged(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1637
MAX_PENALTYPARAM
#define MAX_PENALTYPARAM
Definition: sdpisolver_sdpa.cpp:77
Generated on Tue May 24 2016 11:23:37 for SCIP-SDP by   doxygen 1.8.6