SCIP-SDP  4.0.0
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 programs based on SCIP. */
5 /* */
6 /* Copyright (C) 2011-2013 Discrete Optimization, TU Darmstadt */
7 /* EDOM, FAU Erlangen-Nürnberg */
8 /* 2014-2021 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-2021 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 */
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 
46 #include <assert.h>
47 #ifdef OPENBLAS
48 #include <cblas.h>
49 #endif
50 
51 #include "sdpi/sdpisolver.h"
52 
53 /* turn off warnings for sdpa (doesn't seem to work) */
54 #pragma GCC diagnostic ignored "-Wcast-qual"
55 #pragma GCC diagnostic ignored "-Wpedantic"
56 #include "sdpa_call.h" /* SDPA callable library interface */
57 #pragma GCC diagnostic warning "-Wcast-qual"
58 #pragma GCC diagnostic warning "-Wpedantic"
59 
60 #include "blockmemshell/memory.h" /* for memory allocation */
61 #include "scip/def.h" /* for SCIP_Real, _Bool, ... */
62 #include "scip/pub_misc.h" /* for sorting */
63 #include "sdpi/sdpsolchecker.h" /* to check solution with regards to feasibility tolerance */
64 #include "scip/pub_message.h" /* for debug and error message */
65 
66 /* turn off lint warnings for whole file: */
67 /*lint --e{1784}*/
68 
69 
70 /* local defines */
71 #define GAPTOLCHANGE 1
72 #define FEASTOLCHANGE 1
73 #define PENALTYBOUNDTOL 1E-3
76 #define INFEASFEASTOLCHANGE 0.1
77 #define INFEASMINFEASTOL 1E-9
79 #define MIN_LAMBDASTAR 1e-6
80 #define MAX_LAMBDASTAR 1e8
81 #define LAMBDASTAR_FACTOR 1e0
82 #define LAMBDASTAR_TWOPOINTS TRUE
83 #define LAMBDASTAR_THRESHOLD 1e1
84 #define LAMBDASTAR_LOW 1.5
85 #define LAMBDASTAR_HIGH 1e5
86 #define LAMBDASTAR_DEFAULT 1e2
88 #define MIN_PENALTYPARAM 1e5
89 #define MAX_PENALTYPARAM 1e12
90 #define PENALTYPARAM_FACTOR 1e1
91 #define MAX_MAXPENALTYPARAM 1e15
92 #define MAXPENALTYPARAM_FACTOR 1e6
96 #define BMS_CALL(x) do \
97  { \
98  if( NULL == (x) ) \
99  { \
100  SCIPerrorMessage("No memory in function call.\n"); \
101  return SCIP_NOMEMORY; \
102  } \
103  } \
104  while( FALSE )
105 
107 #define CHECK_IF_SOLVED(sdpisolver) do \
108  { \
109  if (!(sdpisolver->solved)) \
110  { \
111  SCIPerrorMessage("Tried to access solution information for SDP %d ahead of solving!\n", sdpisolver->sdpcounter); \
112  return SCIP_LPERROR; \
113  } \
114  } \
115  while( FALSE )
116 
118 #define CHECK_IF_SOLVED_BOOL(sdpisolver) do \
119  { \
120  if (!(sdpisolver->solved)) \
121  { \
122  SCIPerrorMessage("Tried to access solution information for SDP %d ahead of solving!\n", sdpisolver->sdpcounter); \
123  return FALSE; \
124  } \
125  } \
126  while( FALSE )
127 
128 
130 struct SCIP_SDPiSolver
131 {
132  SCIP_MESSAGEHDLR* messagehdlr;
133  BMS_BLKMEM* blkmem;
134  BMS_BUFMEM* bufmem;
135  SDPA* sdpa;
136  int nvars;
137  int maxnvars;
138  int nactivevars;
139  int* inputtosdpamapper;
142  int* sdpatoinputmapper;
143  SCIP_Real* fixedvarsval;
144  SCIP_Real fixedvarsobjcontr;
145  SCIP_Real* objcoefs;
146  int nvarbounds;
147  int* varboundpos;
149  int* inputtovbmapper;
151  int nsdpalpcons;
152  int nsdpblocks;
153  int maxnsdpblocks;
154  int* maxsdpblocksizes;
155  int maxnlpcons;
156  int* rowmapper;
157  int* rowtoinputmapper;
158  int* inputtoblockmapper;
159  int** blockindmapper;
160  SCIP_Bool solved;
161  SCIP_Bool timelimit;
162  int sdpcounter;
163  int niterations;
164  int nsdpcalls;
165  SCIP_Real epsilon;
166  SCIP_Real gaptol;
167  SCIP_Real feastol;
168  SCIP_Real sdpsolverfeastol;
169  SCIP_Real objlimit;
170  SCIP_Bool sdpinfo;
171  SCIP_Bool penalty;
172  SCIP_Bool feasorig;
174  SCIP_Bool rbound;
175  SCIP_SDPSOLVERSETTING usedsetting;
176  SCIP_Real lambdastar;
177  SCIP_Real* preoptimalsol;
178  SCIP_Real** preoptimalsolx;
179  SCIP_Real* preoptimalsolxlp;
180  SCIP_Bool preoptimalsolexists;
181  SCIP_Real preoptimalgap;
182  int nthreads;
183 };
184 
185 
186 /*
187  * Local Functions
188  */
189 
190 #ifndef NDEBUG
191 
192 static
193 SCIP_Bool isFixed(
194  SCIP_SDPISOLVER* sdpisolver,
195  SCIP_Real lb,
196  SCIP_Real ub
197  )
198 {
199  assert( sdpisolver != NULL );
200  assert( lb < ub + sdpisolver->feastol );
201 
202  return (ub-lb <= sdpisolver->epsilon);
203 }
204 #else
205 #define isFixed(sdpisolver,lb,ub) (ub-lb <= sdpisolver->epsilon)
206 #endif
207 
208 
210 static
211 int calcGrowSize(
212  int initsize,
213  int num
214  )
215 {
216  int oldsize;
217  int size;
218 
219  assert( initsize >= 0 );
220  assert( num >= 0 );
221 
222  /* calculate the size with loop, such that the resulting numbers are always the same (-> block memory) */
223  initsize = MAX(initsize, SCIP_DEFAULT_MEM_ARRAYGROWINIT);
224  size = initsize;
225  oldsize = size - 1;
226 
227  /* second condition checks against overflow */
228  while ( size < num && size > oldsize )
229  {
230  oldsize = size;
231  size = (int)(SCIP_DEFAULT_MEM_ARRAYGROWFAC * size + initsize);
232  }
233 
234  /* if an overflow happened, set the correct value */
235  if ( size <= oldsize )
236  size = num;
237 
238  assert( size >= initsize );
239  assert( size >= num );
240 
241  return size;
242 }
243 
245 static
246 SCIP_RETCODE ensureMappingDataMemory(
247  SCIP_SDPISOLVER* sdpisolver,
248  int nvars,
249  int nsdpblocks,
250  int* sdpblocksizes,
251  int nlpcons,
252  SCIP_Bool usepreoptimalsol
253  )
254 {
255  int newsize;
256  int oldmaxnsdpblocks;
257  int oldmaxnlpcons;
258  int oldmaxnvars;
259  int blocksize;
260  int b;
261 
262  assert( sdpisolver != NULL );
263 
264  oldmaxnlpcons = sdpisolver->maxnlpcons;
265  oldmaxnvars = sdpisolver->maxnvars;
266 
267  /* treat arrays for variables */
268  if ( nvars > sdpisolver->maxnvars )
269  {
270  newsize = calcGrowSize(sdpisolver->maxnvars, nvars);
271 
272  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->inputtosdpamapper, sdpisolver->maxnvars, newsize) );
273  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->sdpatoinputmapper, sdpisolver->maxnvars, newsize) );
274  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->fixedvarsval, sdpisolver->maxnvars, newsize) );
275  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->objcoefs, sdpisolver->maxnvars, newsize) );
276  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->inputtovbmapper, 2 * sdpisolver->maxnvars, 2 * newsize) );
277  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->varboundpos, 2 * sdpisolver->maxnvars, 2 * newsize) );
278 
279  if ( usepreoptimalsol )
280  {
281  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->preoptimalsol, sdpisolver->maxnvars, newsize) );
282  }
283 
284  sdpisolver->maxnvars = newsize;
285  }
286 
287  /* treat arrays for SDP constraints */
288  if ( nsdpblocks > sdpisolver->maxnsdpblocks )
289  {
290  oldmaxnsdpblocks = sdpisolver->maxnsdpblocks;
291  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->inputtoblockmapper, sdpisolver->maxnsdpblocks, nsdpblocks) );
292  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->blockindmapper, sdpisolver->maxnsdpblocks, nsdpblocks) );
293  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->maxsdpblocksizes, sdpisolver->maxnsdpblocks, nsdpblocks) );
294  if ( usepreoptimalsol )
295  {
296  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->preoptimalsolx, sdpisolver->maxnsdpblocks, nsdpblocks) );
297  }
298  sdpisolver->maxnsdpblocks = nsdpblocks;
299 
300  for (b = 0; b < oldmaxnsdpblocks; ++b)
301  {
302  blocksize = sdpblocksizes[b];
303  if ( blocksize > sdpisolver->maxsdpblocksizes[b] )
304  {
305  BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->blockindmapper[b], sdpisolver->maxsdpblocksizes[b], blocksize);
306 
307  if ( usepreoptimalsol )
308  {
309  BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->preoptimalsolx[b], sdpisolver->maxsdpblocksizes[b] * sdpisolver->maxsdpblocksizes[b], blocksize * blocksize);
310  }
311  sdpisolver->maxsdpblocksizes[b] = blocksize;
312  }
313  }
314 
315  for (b = oldmaxnsdpblocks; b < nsdpblocks; ++b)
316  {
317  blocksize = sdpblocksizes[b];
318  BMSallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->blockindmapper[b], blocksize);
319  if ( usepreoptimalsol )
320  {
321  BMSallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->preoptimalsolx[b], blocksize * blocksize);
322  }
323  sdpisolver->maxsdpblocksizes[b] = blocksize;
324  }
325  }
326 
327  /* treat arrays for LP constraints */
328  if ( nlpcons > sdpisolver->maxnlpcons )
329  {
330  newsize = calcGrowSize(sdpisolver->maxnlpcons, nlpcons);
331  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->rowmapper, 2 * sdpisolver->maxnlpcons, 2 * newsize) ); /*lint !e647*/
332  BMS_CALL( BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->rowtoinputmapper, 2 * sdpisolver->maxnlpcons, 2 * newsize) );
333  sdpisolver->maxnlpcons = newsize;
334  }
335 
336  /* special case: treat preoptimalsolxlp */
337  if ( usepreoptimalsol && ( nlpcons > oldmaxnlpcons || nvars > oldmaxnvars ) )
338  {
339  assert( nlpcons <= sdpisolver->maxnlpcons );
340  assert( nvars <= sdpisolver->maxnvars );
341 
342  BMSreallocBlockMemoryArray(sdpisolver->blkmem, &sdpisolver->preoptimalsolxlp, 2 * (oldmaxnlpcons + oldmaxnvars), 2 * (sdpisolver->maxnlpcons + sdpisolver->maxnvars));
343  }
344 
345  return SCIP_OKAY;
346 }
347 
348 
352 static
353 SCIP_RETCODE checkFeastolAndResolve(
354  SCIP_SDPISOLVER* sdpisolver,
355  SCIP_Real penaltyparam,
356  int nvars,
357  SCIP_Real* lb,
358  SCIP_Real* ub,
359  int nsdpblocks,
360  int* sdpblocksizes,
361  int* sdpnblockvars,
362  int sdpconstnnonz,
363  int* sdpconstnblocknonz,
365  int** sdpconstrow,
366  int** sdpconstcol,
367  SCIP_Real** sdpconstval,
368  int sdpnnonz,
369  int** sdpnblockvarnonz,
371  int** sdpvar,
373  int*** sdprow,
374  int*** sdpcol,
375  SCIP_Real*** sdpval,
376  int** indchanges,
378  int* nremovedinds,
379  int* blockindchanges,
380  int nlpcons,
381  SCIP_Real* lplhs,
382  SCIP_Real* lprhs,
383  int lpnnonz,
384  int* lprow,
385  int* lpcol,
386  SCIP_Real* lpval,
387  SCIP_Real* feastol,
388  SCIP_Real* gaptol
389  )
390 {
391 #ifdef SCIP_DEBUG
392  char phase_string[15];
393 #endif
394 
395  assert( feastol != NULL );
396  assert( gaptol != NULL );
397 
398  if ( penaltyparam >= sdpisolver->epsilon )
399  return SCIP_OKAY;
400 
401  while ( SCIPsdpiSolverIsAcceptable(sdpisolver) && SCIPsdpiSolverIsDualFeasible(sdpisolver) && *feastol >= INFEASMINFEASTOL && *gaptol >= INFEASMINFEASTOL )
402  {
403  SCIP_Real* solvector;
404  SCIP_Real primalobj;
405  SCIP_Real dualobj;
406  SCIP_Bool infeasible;
407  SCIP_Bool solveagain = FALSE;
408  int nvarspointer;
409 
410  /* get current solution */
411  BMS_CALL( BMSallocBufferMemoryArray(sdpisolver->bufmem, &solvector, nvars) ); /*lint !e530*/
412  nvarspointer = nvars;
413  SCIP_CALL( SCIPsdpiSolverGetSol(sdpisolver, NULL, solvector, &nvarspointer) );
414  assert( nvarspointer == nvars );
415 
416  /* check the solution for feasibility with regards to our tolerance */
417  SCIP_CALL( SCIPsdpSolcheckerCheck(sdpisolver->bufmem, nvars, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
418  sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval,
419  indchanges, nremovedinds, blockindchanges, nlpcons, lplhs, lprhs, lpnnonz, lprow, lpcol, lpval,
420  solvector, sdpisolver->feastol, sdpisolver->epsilon, &infeasible) );
421  BMSfreeBufferMemoryArray(sdpisolver->bufmem, &solvector);
422 
423  if ( infeasible )
424  {
425  *feastol *= INFEASFEASTOLCHANGE;
426  if ( *feastol >= INFEASMINFEASTOL )
427  {
428  SCIPdebugMessage("Solution feasible for SDPA but outside feasibility tolerance, changing SDPA feasibility tolerance to %g.\n", *feastol);
429  sdpisolver->sdpa->setParameterEpsilonDash(*feastol);
430  solveagain = TRUE;
431  }
432  }
433 
434  /* check whether duality gap is small enough */
435  dualobj = sdpisolver->sdpa->getPrimalObj();
436  primalobj = sdpisolver->sdpa->getDualObj();
437  if ( REALABS(dualobj - primalobj) >= sdpisolver->gaptol )
438  {
439  infeasible = TRUE;
440  *gaptol *= INFEASFEASTOLCHANGE;
441  if ( *gaptol >= INFEASMINFEASTOL )
442  {
443  SCIPdebugMessage("Solution feasible, but duality gap is too large, changing SDPA gap tolerance to %g.\n", *gaptol);
444  sdpisolver->sdpa->setParameterEpsilonStar(*gaptol);
445  solveagain = TRUE;
446  }
447  }
448 
449  if ( solveagain )
450  {
451 #ifdef SCIP_MORE_DEBUG
452  sdpisolver->sdpa->printParameters(stdout);
453 #endif
454  sdpisolver->sdpa->setInitPoint(false);
455 #ifdef SDPA_RESETPARAMS
456  sdpisolver->sdpa->resetParameters();
457 #else
458  sdpisolver->sdpa->initializeSolve();
459 #endif
460  sdpisolver->sdpa->solve();
461 
462  /* update number of SDP-iterations and -calls */
463  sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
464  ++sdpisolver->nsdpcalls;
465 
466 #ifdef SCIP_DEBUG
467  /* print the phase value , i.e. whether solving was successfull */
468  sdpisolver->sdpa->getPhaseString((char*)phase_string);
469  SCIPdebugMessage("SDPA solving finished with status %s (primal and dual here are switched in contrast to our formulation)\n", phase_string);
470 #endif
471  }
472  else
473  {
474  if ( infeasible )
475  {
476  sdpisolver->solved = FALSE;
477  SCIPmessagePrintInfo(sdpisolver->messagehdlr, "SDPA failed to reach required feasibility tolerance!\n");
478  }
479  break;
480  }
481  }
482 
483  return SCIP_OKAY;
484 }
485 
486 /*
487  * Miscellaneous Methods
488  */
489 
495 const char* SCIPsdpiSolverGetSolverName(
496  void
497  )
498 {/*lint !e1784*/
499  return "SDPA"; /* version number can only be printed to file/stdout */
500 }
501 
503 const char* SCIPsdpiSolverGetSolverDesc(
504  void
505  )
506 {/*lint !e1784*/
507  return "Primal-dual Interior Point Solver for SDPs developed by K. Fujisawa et al. (sdpa.sourceforge.net)";
508 }
509 
516 void* SCIPsdpiSolverGetSolverPointer(
517  SCIP_SDPISOLVER* sdpisolver
518  )
519 {/*lint !e1784*/
520  assert( sdpisolver != NULL );
521  return (void*) sdpisolver->sdpa;
522 }
523 
525 int SCIPsdpiSolverGetDefaultSdpiSolverNpenaltyIncreases(
526  void
527  )
528 {
529  return 2;
530 }
531 
533 SCIP_Bool SCIPsdpiSolverDoesWarmstartNeedPrimal(
534  void
535  )
536 {
537  return TRUE;
538 }
539 
543 /*
544  * SDPI Creation and Destruction Methods
545  */
546 
551 SCIP_RETCODE SCIPsdpiSolverCreate(
552  SCIP_SDPISOLVER** sdpisolver,
553  SCIP_MESSAGEHDLR* messagehdlr,
554  BMS_BLKMEM* blkmem,
555  BMS_BUFMEM* bufmem
556  )
557 {/*lint !e1784*/
558  assert( sdpisolver != NULL );
559  assert( blkmem != NULL );
560  assert( bufmem != NULL );
561 
562  SCIPdebugMessage("Calling SCIPsdpiCreate \n");
563 
564  BMS_CALL( BMSallocBlockMemory(blkmem, sdpisolver) );
565 
566  (*sdpisolver)->messagehdlr = messagehdlr;
567  (*sdpisolver)->blkmem = blkmem;
568  (*sdpisolver)->bufmem = bufmem;
569 
570  /* this will be properly initialized then calling solve */
571  (*sdpisolver)->sdpa = NULL;
572 
573  (*sdpisolver)->nvars = 0;
574  (*sdpisolver)->maxnvars = 0;
575  (*sdpisolver)->nactivevars = 0;
576  (*sdpisolver)->nsdpblocks = 0;
577  (*sdpisolver)->maxnsdpblocks = 0;
578  (*sdpisolver)->inputtosdpamapper = NULL;
579  (*sdpisolver)->sdpatoinputmapper = NULL;
580  (*sdpisolver)->inputtoblockmapper = NULL;
581  (*sdpisolver)->blockindmapper = NULL;
582  (*sdpisolver)->maxsdpblocksizes = NULL;
583  (*sdpisolver)->maxnlpcons = 0;
584  (*sdpisolver)->rowmapper = NULL;
585  (*sdpisolver)->rowtoinputmapper = NULL;
586  (*sdpisolver)->fixedvarsval = NULL;
587  (*sdpisolver)->fixedvarsobjcontr = 0.0;
588  (*sdpisolver)->objcoefs = NULL;
589  (*sdpisolver)->nvarbounds = 0;
590  (*sdpisolver)->varboundpos = NULL;
591  (*sdpisolver)->inputtovbmapper = NULL;
592  (*sdpisolver)->solved = FALSE;
593  (*sdpisolver)->timelimit = FALSE;
594  (*sdpisolver)->sdpcounter = 0;
595  (*sdpisolver)->niterations = 0;
596  (*sdpisolver)->nsdpcalls = 0;
597  (*sdpisolver)->nsdpalpcons = 0;
598 
599  (*sdpisolver)->epsilon = 1e-9;
600  (*sdpisolver)->gaptol = 1e-6;
601  (*sdpisolver)->feastol = 1e-6;
602  (*sdpisolver)->sdpsolverfeastol = 1e-6;
603  (*sdpisolver)->objlimit = SCIPsdpiSolverInfinity(*sdpisolver);
604  (*sdpisolver)->sdpinfo = FALSE;
605  (*sdpisolver)->usedsetting = SCIP_SDPSOLVERSETTING_UNSOLVED;
606  (*sdpisolver)->lambdastar = LAMBDASTAR_DEFAULT;
607 
608  (*sdpisolver)->preoptimalsol = NULL;
609  (*sdpisolver)->preoptimalsolx = NULL;
610  (*sdpisolver)->preoptimalsolxlp = NULL;
611  (*sdpisolver)->preoptimalsolexists = FALSE;
612  (*sdpisolver)->preoptimalgap = -1.0;
613 
614  (*sdpisolver)->nthreads = -1;
615 
616  return SCIP_OKAY;
617 }
618 
620 SCIP_RETCODE SCIPsdpiSolverFree(
621  SCIP_SDPISOLVER** sdpisolver
622  )
623 {/*lint !e1784*/
624  int b;
625  int blocksize;
626 
627  assert( sdpisolver != NULL );
628  assert( *sdpisolver != NULL );
629 
630  SCIPdebugMessage("Freeing SDPISolver\n");
631 
632  for (b = 0; b < (*sdpisolver)->maxnsdpblocks; b++)
633  {
634  blocksize = (*sdpisolver)->maxsdpblocksizes[b];
635  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &((*sdpisolver)->blockindmapper[b]), blocksize);
636 
637  if ( (*sdpisolver)->preoptimalsolx != NULL )
638  {
639  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &((*sdpisolver)->preoptimalsolx[b]), blocksize * blocksize);
640  }
641  }
642  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->preoptimalsolxlp, 2 * ((*sdpisolver)->maxnlpcons + (*sdpisolver)->maxnvars));
643 
644  /* free LP block of preoptimal solution X array */
645  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->preoptimalsolx, (*sdpisolver)->maxnsdpblocks);
646  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->blockindmapper, (*sdpisolver)->maxnsdpblocks);
647  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->maxsdpblocksizes, (*sdpisolver)->maxnsdpblocks);
648 
649  /* free SDPA object using destructor and free memory via blockmemshell */
650  if ( (*sdpisolver)->sdpa != NULL)
651  delete (*sdpisolver)->sdpa;
652 
653  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->preoptimalsol, (*sdpisolver)->maxnvars);
654  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->inputtoblockmapper, (*sdpisolver)->maxnsdpblocks);
655 
656  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->rowtoinputmapper, 2 * (*sdpisolver)->maxnlpcons);
657  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->rowmapper, 2 * (*sdpisolver)->maxnlpcons);
658 
659  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->varboundpos, 2 * (*sdpisolver)->maxnvars);
660  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->inputtovbmapper, 2 * (*sdpisolver)->maxnvars);
661  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->inputtosdpamapper, (*sdpisolver)->maxnvars);
662  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->sdpatoinputmapper, (*sdpisolver)->maxnvars);
663  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->fixedvarsval, (*sdpisolver)->maxnvars);
664  BMSfreeBlockMemoryArrayNull((*sdpisolver)->blkmem, &(*sdpisolver)->objcoefs, (*sdpisolver)->maxnvars);
665  BMSfreeBlockMemory((*sdpisolver)->blkmem, sdpisolver);
666 
667  return SCIP_OKAY;
668 }
669 
671 SCIP_RETCODE SCIPsdpiSolverIncreaseCounter(
672  SCIP_SDPISOLVER* sdpisolver
673  )
674 {/*lint !e1784*/
675  assert( sdpisolver != NULL );
676 
677  sdpisolver->sdpcounter++;
678 
679  return SCIP_OKAY;
680 }
681 
683 SCIP_RETCODE SCIPsdpiSolverResetCounter(
684  SCIP_SDPISOLVER* sdpisolver
685  )
686 {/*lint !e1784*/
687  assert( sdpisolver != NULL );
688 
689  SCIPdebugMessage("Resetting counter of SDP-Interface from %d to 0.\n", sdpisolver->sdpcounter);
690  sdpisolver->sdpcounter = 0;
691 
692  return SCIP_OKAY;
693 }
694 
698 /*
699  * Solving Methods
700  */
701 
720 SCIP_RETCODE SCIPsdpiSolverLoadAndSolve(
721  SCIP_SDPISOLVER* sdpisolver,
722  int nvars,
723  SCIP_Real* obj,
724  SCIP_Real* lb,
725  SCIP_Real* ub,
726  int nsdpblocks,
727  int* sdpblocksizes,
728  int* sdpnblockvars,
729  int sdpconstnnonz,
730  int* sdpconstnblocknonz,
732  int** sdpconstrow,
733  int** sdpconstcol,
734  SCIP_Real** sdpconstval,
735  int sdpnnonz,
736  int** sdpnblockvarnonz,
738  int** sdpvar,
740  int*** sdprow,
741  int*** sdpcol,
742  SCIP_Real*** sdpval,
743  int** indchanges,
745  int* nremovedinds,
746  int* blockindchanges,
747  int nremovedblocks,
748  int nlpcons,
749  SCIP_Real* lplhs,
750  SCIP_Real* lprhs,
751  int lpnnonz,
752  int* lprow,
753  int* lpcol,
754  SCIP_Real* lpval,
755  SCIP_Real* starty,
756  int* startZnblocknonz,
758  int** startZrow,
760  int** startZcol,
762  SCIP_Real** startZval,
764  int* startXnblocknonz,
766  int** startXrow,
768  int** startXcol,
770  SCIP_Real** startXval,
772  SCIP_SDPSOLVERSETTING startsettings,
774  SCIP_Real timelimit,
775  SDPI_CLOCK* usedsdpitime
776  )
777 {/*lint !e1784*/
778  return SCIPsdpiSolverLoadAndSolveWithPenalty(sdpisolver, 0.0, TRUE, FALSE, nvars, obj, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
779  sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
780  nremovedinds, blockindchanges, nremovedblocks, nlpcons, lplhs, lprhs, lpnnonz, lprow, lpcol, lpval,
781  starty, startZnblocknonz, startZrow, startZcol, startZval, startXnblocknonz, startXrow, startXcol, startXval, startsettings,
782  timelimit, usedsdpitime, NULL, NULL);
783 }
784 
807 SCIP_RETCODE SCIPsdpiSolverLoadAndSolveWithPenalty(
808  SCIP_SDPISOLVER* sdpisolver,
809  SCIP_Real penaltyparam,
810  SCIP_Bool withobj,
811  SCIP_Bool rbound,
812  int nvars,
813  SCIP_Real* obj,
814  SCIP_Real* lb,
815  SCIP_Real* ub,
816  int nsdpblocks,
817  int* sdpblocksizes,
818  int* sdpnblockvars,
819  int sdpconstnnonz,
820  int* sdpconstnblocknonz,
822  int** sdpconstrow,
823  int** sdpconstcol,
824  SCIP_Real** sdpconstval,
825  int sdpnnonz,
826  int** sdpnblockvarnonz,
828  int** sdpvar,
830  int*** sdprow,
831  int*** sdpcol,
832  SCIP_Real*** sdpval,
833  int** indchanges,
835  int* nremovedinds,
836  int* blockindchanges,
837  int nremovedblocks,
838  int nlpcons,
839  SCIP_Real* lplhs,
840  SCIP_Real* lprhs,
841  int lpnnonz,
842  int* lprow,
843  int* lpcol,
844  SCIP_Real* lpval,
845  SCIP_Real* starty,
846  int* startZnblocknonz,
848  int** startZrow,
850  int** startZcol,
852  SCIP_Real** startZval,
854  int* startXnblocknonz,
856  int** startXrow,
858  int** startXcol,
860  SCIP_Real** startXval,
862  SCIP_SDPSOLVERSETTING startsettings,
864  SCIP_Real timelimit,
865  SDPI_CLOCK* usedsdpitime,
866  SCIP_Bool* feasorig,
868  SCIP_Bool* penaltybound
870  )
871 {/*lint !e1784*/
872  SCIP_Real feastol;
873  SCIP_Real gaptol;
874  SCIP_Real solvertimelimit;
875  SCIP_Bool usepreoptimalsol;
876  int nfixedvars;
877  bool checkinput; /* should the input be checked for consistency in SDPA? */
878  int nsdpablocks;
879  int nlpineqs;
880  int ind;
881  int i;
882  int k;
883  int b;
884 
885 #ifdef SCIP_DEBUG
886  char phase_string[15];
887 #endif
888 
889  assert( sdpisolver != NULL );
890  assert( penaltyparam > -1 * sdpisolver->epsilon );
891  assert( penaltyparam < sdpisolver->epsilon || ( feasorig != NULL ) );
892  assert( nvars > 0 );
893  assert( obj != NULL );
894  assert( lb != NULL );
895  assert( ub != NULL );
896  assert( nsdpblocks >= 0 );
897  assert( nsdpblocks == 0 || sdpblocksizes != NULL );
898  assert( nsdpblocks == 0 || sdpnblockvars != NULL );
899  assert( sdpconstnnonz >= 0 );
900  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstnblocknonz != NULL );
901  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstrow != NULL );
902  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstcol != NULL );
903  assert( nsdpblocks == 0 || sdpconstnnonz == 0 || sdpconstval != NULL );
904  assert( sdpnnonz >= 0 );
905  assert( nsdpblocks == 0 || sdpnblockvarnonz != NULL );
906  assert( nsdpblocks == 0 || sdpvar != NULL );
907  assert( nsdpblocks == 0 || sdprow != NULL );
908  assert( nsdpblocks == 0 || sdpcol != NULL );
909  assert( nsdpblocks == 0 || sdpval != NULL );
910  assert( nsdpblocks == 0 || indchanges != NULL );
911  assert( nsdpblocks == 0 || nremovedinds != NULL );
912  assert( nsdpblocks == 0 || blockindchanges != NULL );
913  assert( 0 <= nremovedblocks && nremovedblocks <= nsdpblocks );
914  assert( nlpcons >= 0 );
915  assert( nlpcons == 0 || lplhs != NULL );
916  assert( nlpcons == 0 || lprhs != NULL );
917  assert( lpnnonz >= 0 );
918  assert( nlpcons == 0 || lprow != NULL );
919  assert( nlpcons == 0 || lpcol != NULL );
920  assert( nlpcons == 0 || lpval != NULL );
921  assert( startZnblocknonz == NULL || startZrow != NULL );
922  assert( startZnblocknonz == NULL || startZcol != NULL );
923  assert( startZnblocknonz == NULL || startZval != NULL );
924  assert( startXnblocknonz == NULL || startXrow != NULL );
925  assert( startXnblocknonz == NULL || startXcol != NULL );
926  assert( startXnblocknonz == NULL || startXval != NULL );
927 
928  sdpisolver->niterations = 0;
929  sdpisolver->nsdpcalls = 0;
930  sdpisolver->feasorig = FALSE;
931 
932  /* compute the time limit to set for the solver */
933  solvertimelimit = timelimit;
934  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, solvertimelimit) )
935  solvertimelimit -= SDPIclockGetTime(usedsdpitime);
936 
937  if ( solvertimelimit <= 0.0 )
938  {
939  sdpisolver->solved = FALSE;
940  sdpisolver->timelimit = TRUE;
941  return SCIP_OKAY;
942  }
943  else
944  sdpisolver->timelimit = FALSE;
945 
946 #ifndef NDEBUG
947  checkinput = true;
948 #else
949  checkinput = false;
950 #endif
951 
952  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_UNSOLVED;
953 
954  /* 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
955  * same SDP) */
956  if ( penaltyparam < sdpisolver->epsilon )
957  SCIPdebugMessage("Inserting Data into SDPA for SDP (%d) \n", ++sdpisolver->sdpcounter);
958  else
959  SCIPdebugMessage("Inserting Data again into SDPA for SDP (%d) \n", sdpisolver->sdpcounter);
960 
961  /* free old SDPA solver */
962  if ( sdpisolver->sdpa != 0 )
963  {
964  /* if the SDPA solver has already been created, clear the current problem instance */
965  delete sdpisolver->sdpa;
966  }
967  sdpisolver->sdpa = new SDPA();
968  assert( sdpisolver->sdpa != 0 );
969 
970  /* set number of threads (do this early, since this might affect factorization) */
971  if ( sdpisolver->nthreads > 0 )
972  {
973 #ifdef OPENBLAS
974  openblas_set_num_threads(sdpisolver->nthreads);
975 #endif
976  sdpisolver->sdpa->setNumThreads(sdpisolver->nthreads);
977  }
978 
979  /* set the penalty and rbound flags accordingly */
980  sdpisolver->penalty = (penaltyparam < sdpisolver->epsilon) ? FALSE : TRUE;
981  sdpisolver->rbound = rbound;
982 
983  if ( sdpisolver->preoptimalgap >= 0.0 )
984  usepreoptimalsol = TRUE;
985  else
986  usepreoptimalsol = FALSE;
987 
988  /* ensure memory for working space */
989  SCIP_CALL( ensureMappingDataMemory(sdpisolver, nvars, nsdpblocks, sdpblocksizes, nlpcons, usepreoptimalsol) );
990 
991  assert( sdpisolver->nvars <= sdpisolver->maxnvars );
992  sdpisolver->nvars = nvars;
993  sdpisolver->nactivevars = 0;
994  sdpisolver->nvarbounds = 0;
995  sdpisolver->preoptimalsolexists = FALSE;
996  sdpisolver->nsdpblocks = nsdpblocks;
997  sdpisolver->nsdpalpcons = nlpcons;
998  nsdpablocks = nsdpblocks - nremovedblocks;
999  nfixedvars = 0;
1000 
1001  /* find fixed variables and set up mapping functions */
1002  sdpisolver->fixedvarsobjcontr = 0.0;
1003  for (i = 0; i < nvars; i++)
1004  {
1005  if ( isFixed(sdpisolver, lb[i], ub[i]) )
1006  {
1007  sdpisolver->fixedvarsobjcontr += obj[i] * lb[i]; /* this is the value this fixed variable contributes to the objective */
1008  sdpisolver->fixedvarsval[nfixedvars] = lb[i]; /* if lb=ub, then this is the value the variable will have in every solution */
1009  nfixedvars++;
1010  sdpisolver->inputtosdpamapper[i] = -nfixedvars;
1011  SCIPdebugMessage("Fixing variable %d locally to %g in SDPA.\n", i, lb[i]);
1012  sdpisolver->inputtovbmapper[2 * i] = -1;
1013  sdpisolver->inputtovbmapper[2 * i + 1] = -1;
1014  }
1015  else
1016  {
1017 #ifdef SCIP_MORE_DEBUG
1018  SCIPdebugMessage("Variable %d becomes variable %d in SDPA.\n", i, sdpisolver->nactivevars + 1);
1019 #endif
1020  sdpisolver->sdpatoinputmapper[sdpisolver->nactivevars] = i;
1021  sdpisolver->inputtosdpamapper[i] = sdpisolver->nactivevars + 1; /* sdpa starts counting at 1 */
1022  sdpisolver->objcoefs[sdpisolver->nactivevars] = obj[i];
1023 
1024  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, lb[i]) )
1025  {
1026  sdpisolver->varboundpos[sdpisolver->nvarbounds] = -(sdpisolver->nactivevars + 1); /* negative sign means lower bound, +1 because of SDPA */
1027  sdpisolver->inputtovbmapper[2 * i] = sdpisolver->nvarbounds++;
1028  }
1029 
1030  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, ub[i]) )
1031  {
1032  sdpisolver->varboundpos[sdpisolver->nvarbounds] = +(sdpisolver->nactivevars + 1); /* positive sign means upper bound, +1 because of SDPA */
1033  sdpisolver->inputtovbmapper[2 * i + 1] = sdpisolver->nvarbounds++;
1034  }
1035 
1036  sdpisolver->nactivevars++;
1037  }
1038  }
1039  assert( sdpisolver->nactivevars + nfixedvars == sdpisolver->nvars );
1040  assert( sdpisolver->nvarbounds <= 2 * sdpisolver->nactivevars );
1041 
1042  /* if we use a penalty formulation, we need the constraint r >= 0 */
1043  if ( penaltyparam >= sdpisolver->epsilon && rbound )
1044  sdpisolver->nvarbounds++;
1045 
1046  /* if we want to solve without objective, we reset fixedvarsobjcontr */
1047  if ( ! withobj )
1048  sdpisolver->fixedvarsobjcontr = 0.0;
1049 
1050  /* determine total number of sides in LP-constraints */
1051  nlpineqs = 0;
1052  for (i = 0; i < nlpcons; i++)
1053  {
1054  assert( lplhs[i] < SCIPsdpiSolverInfinity(sdpisolver) );
1055  if ( lplhs[i] > - SCIPsdpiSolverInfinity(sdpisolver) )
1056  {
1057  sdpisolver->rowmapper[2 * i] = nlpineqs + 1;
1058  sdpisolver->rowtoinputmapper[nlpineqs] = 2 * i;
1059  ++nlpineqs;
1060  }
1061  else
1062  sdpisolver->rowmapper[2 * i] = -1;
1063 
1064  assert( lprhs[i] > -SCIPsdpiSolverInfinity(sdpisolver) );
1065  if ( lprhs[i] < SCIPsdpiSolverInfinity(sdpisolver) )
1066  {
1067  sdpisolver->rowmapper[2*i + 1] = nlpineqs + 1;
1068  sdpisolver->rowtoinputmapper[nlpineqs] = 2 * i + 1;
1069  ++nlpineqs;
1070  }
1071  else
1072  sdpisolver->rowmapper[2*i + 1] = -1;
1073  }
1074  assert( nlpineqs <= 2*nlpcons );
1075 
1076  /* compute blockmapper and blockindmapper */
1077  ind = 0;
1078  assert( nsdpblocks <= sdpisolver->maxnsdpblocks );
1079  for (b = 0; b < nsdpblocks; ++b)
1080  {
1081  if ( blockindchanges[b] >= 0 )
1082  {
1083  int blockind = 0;
1084 
1085  sdpisolver->inputtoblockmapper[b] = ind + 1; /* +1 because SDPA start counting at one */
1086 
1087  for (i = 0; i < sdpblocksizes[b]; i++)
1088  {
1089  if ( indchanges[b][i] >= 0 )
1090  {
1091  assert( ind <= nsdpblocks );
1092  sdpisolver->blockindmapper[ind][blockind] = i;
1093  assert( i - indchanges[b][i] == blockind );
1094  blockind++;
1095  }
1096  }
1097  assert( blockind == sdpblocksizes[b] - nremovedinds[b] );
1098  ind++;
1099  }
1100  else
1101  sdpisolver->inputtoblockmapper[b] = -1;
1102  }
1103 
1104  /* initialize blockstruct */
1105  if ( penaltyparam < sdpisolver->epsilon ) /* we initialize this with an exact 0.0 in solve without penalty */
1106  sdpisolver->sdpa->inputConstraintNumber((long long) sdpisolver->nactivevars);
1107  else
1108  sdpisolver->sdpa->inputConstraintNumber((long long) sdpisolver->nactivevars + 1); /* the additional variable is r for the penalty formulation */
1109 
1110  /* if there are any LP-cons/variable-bounds, we get an extra block for those */
1111  if ( nlpineqs + sdpisolver->nvarbounds > 0 )
1112  sdpisolver->sdpa->inputBlockNumber((long long) (nsdpablocks + 1));
1113  else
1114  sdpisolver->sdpa->inputBlockNumber((long long) (nsdpablocks));
1115 
1116  /* initialize SDP blocks */
1117  for (b = 0; b < nsdpblocks; ++b)
1118  {
1119  if ( blockindchanges[b] >= 0 )
1120  {
1121  /* +1 because of SDPA counting */
1122  SCIPdebugMessage("adding block %d to SDPA as block %d with size %d.\n", b, b - blockindchanges[b] + 1, sdpblocksizes[b] - nremovedinds[b]);
1123  sdpisolver->sdpa->inputBlockSize((long long) (b - blockindchanges[b] + 1), (long long) (sdpblocksizes[b] - nremovedinds[b]));
1124  sdpisolver->sdpa->inputBlockType((long long) (b - blockindchanges[b] + 1), SDPA::SDP);
1125  }
1126  }
1127 
1128  /* initialize LP block */
1129  if ( nlpineqs + sdpisolver->nvarbounds > 0 )
1130  {
1131  /* the last block is the LP block, the size has a negative sign */
1132  SCIPdebugMessage("adding LP block to SDPA as block %d with size %d\n", nsdpablocks + 1, nlpineqs + sdpisolver->nvarbounds);
1133  sdpisolver->sdpa->inputBlockSize((long long) (nsdpablocks + 1), (long long) -(nlpineqs + sdpisolver->nvarbounds));
1134  sdpisolver->sdpa->inputBlockType((long long) (nsdpablocks + 1), SDPA::LP);
1135  }
1136 
1137  /* initialize space */
1138  sdpisolver->sdpa->initializeUpperTriangleSpace();
1139 
1140  /* set objective values */
1141  if ( withobj )
1142  {
1143  for (i = 0; i < sdpisolver->nactivevars; i++)
1144  {
1145  /* insert objective value, SDPA counts from 1 to n instead of 0 to n-1 */
1146  sdpisolver->sdpa->inputCVec((long long) (i + 1), obj[sdpisolver->sdpatoinputmapper[i]]);
1147  }
1148  }
1149 
1150  if ( penaltyparam >= sdpisolver->epsilon )
1151  sdpisolver->sdpa->inputCVec((long long) (sdpisolver->nactivevars + 1), penaltyparam); /* set the objective of the additional var to penaltyparam */
1152 
1153  /* inserting non-constant SDP constraint matrices */
1154  if ( sdpnnonz > 0 )
1155  {
1156  int blockvar;
1157  int v;
1158 
1159  assert( nsdpblocks > 0 );
1160  assert( sdpnblockvarnonz != NULL );
1161  assert( sdpnblockvars != NULL );
1162  assert( sdpcol != NULL );
1163  assert( sdprow != NULL );
1164  assert( sdpval != NULL );
1165  assert( sdpvar != NULL );
1166  assert( indchanges != NULL );
1167  assert( nremovedinds != NULL );
1168 
1169  for (b = 0; b < nsdpblocks; ++b)
1170  {
1171  /* if the block has no entries, we skip it */
1172  if ( blockindchanges[b] < 0 )
1173  continue;
1174 
1175  assert( b - blockindchanges[b] + 1 <= nsdpablocks );
1176 
1177  /* iterate over all variables in this block */
1178  for (blockvar = 0; blockvar < sdpnblockvars[b]; blockvar++)
1179  {
1180  v = sdpisolver->inputtosdpamapper[sdpvar[b][blockvar]];
1181 
1182  /* check if the variable is active */
1183  if ( v > 0 )
1184  {
1185  assert( v <= sdpisolver->nactivevars );
1186 
1187  for (k = 0; k < sdpnblockvarnonz[b][blockvar]; k++)
1188  {
1189  /* rows and cols with active nonzeros should not be removed */
1190  assert( 0 <= indchanges[b][sdprow[b][blockvar][k]] && indchanges[b][sdprow[b][blockvar][k]] <= sdprow[b][blockvar][k] );
1191  assert( 0 <= indchanges[b][sdpcol[b][blockvar][k]] && indchanges[b][sdpcol[b][blockvar][k]] <= sdpcol[b][blockvar][k] );
1192 
1193  assert( 0 <= sdprow[b][blockvar][k] && sdprow[b][blockvar][k] < sdpblocksizes[b] );
1194  assert( 0 <= sdpcol[b][blockvar][k] && sdpcol[b][blockvar][k] < sdpblocksizes[b] );
1195 
1196  /* rows and columns start with 1 in SDPA, so we have to add 1 to the indices */
1197  sdpisolver->sdpa->inputElement((long long) v, (long long) (b - blockindchanges[b] + 1),
1198  (long long) (sdpcol[b][blockvar][k] - indchanges[b][sdpcol[b][blockvar][k]] + 1),
1199  (long long) (sdprow[b][blockvar][k] - indchanges[b][sdprow[b][blockvar][k]] + 1),
1200  sdpval[b][blockvar][k], checkinput);
1201  }
1202  }
1203  }
1204 
1205  /* insert the identity matrix if we are using a penalty formulation */
1206  if ( penaltyparam >= sdpisolver->epsilon )
1207  {
1208  for (i = 0; i < sdpblocksizes[b] - nremovedinds[b]; i++)
1209  {
1210  sdpisolver->sdpa->inputElement((long long) (sdpisolver->nactivevars + 1), (long long) (b - blockindchanges[b] + 1),
1211  (long long) (i + 1), (long long) (i + 1), 1.0, checkinput);
1212  }
1213  }
1214  }
1215  }
1216 
1217  /* inserting constant matrix */
1218  if ( sdpconstnnonz > 0 )
1219  {
1220  assert( nsdpblocks > 0 );
1221  assert( sdpconstnblocknonz!= NULL );
1222  assert( sdpconstcol != NULL );
1223  assert( sdpconstrow != NULL );
1224  assert( sdpconstval != NULL );
1225 
1226  for (b = 0; b < nsdpblocks; ++b)
1227  {
1228  if ( blockindchanges[b] < 0 )
1229  continue;
1230 
1231  for (k = 0; k < sdpconstnblocknonz[b]; k++)
1232  {
1233  /* rows and cols with active nonzeros should not be removed */
1234  assert( -1 < indchanges[b][sdpconstrow[b][k]] && indchanges[b][sdpconstrow[b][k]] <= sdpconstrow[b][k] );
1235  assert( -1 < indchanges[b][sdpconstcol[b][k]] && indchanges[b][sdpconstcol[b][k]] <= sdpconstcol[b][k] );
1236 
1237  assert( 0 <= sdpconstrow[b][k] && sdpconstrow[b][k] < sdpblocksizes[b] );
1238  assert( 0 <= sdpconstcol[b][k] && sdpconstcol[b][k] < sdpblocksizes[b] );
1239 
1240  /* rows and columns start with 1 in SDPA, so we have to add 1 to the indices, the constant matrix is given as variable 0 */
1241  sdpisolver->sdpa->inputElement((long long) 0, (long long) (b - blockindchanges[b] + 1),
1242  (long long) (sdpconstcol[b][k] - indchanges[b][sdpconstcol[b][k]] + 1),
1243  (long long) (sdpconstrow[b][k] - indchanges[b][sdpconstrow[b][k]] + 1),
1244  sdpconstval[b][k], checkinput);
1245  }
1246  }
1247  }
1248 
1249  /* inserting LP nonzeros */
1250  if ( lpnnonz > 0 )
1251  {
1252  SCIP_Bool lhspresent = FALSE;
1253  SCIP_Bool rhspresent = FALSE;
1254  long long lhsrhscnt = 1; /* 1 because SDPA starts with 1 */
1255  int lastrow = -1;
1256  int v;
1257 
1258  for (i = 0; i < lpnnonz; i++)
1259  {
1260  assert( i == 0 || lprow[i-1] <= lprow[i] ); /* rows should be sorted */
1261  assert( 0 <= lprow[i] && lprow[i] < nlpcons );
1262  assert( 0 <= lpcol[i] && lpcol[i] < nvars );
1263  assert( REALABS(lpval[i]) > sdpisolver->epsilon );
1264 
1265  /* if we found a new row */
1266  if ( lprow[i] > lastrow )
1267  {
1268  /* move counter forward */
1269  if ( lhspresent )
1270  {
1271  assert( 0 <= lastrow && lastrow < nlpcons );
1272 
1273  /* LP constraints are added as diagonal entries of the last block, left-hand-side is added as variable zero */
1274  sdpisolver->sdpa->inputElement(0LL, (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, lplhs[lastrow], checkinput);
1275 
1276  /* if we use a penalty formulation, add the r * Identity entry */
1277  if ( penaltyparam >= sdpisolver->epsilon )
1278  {
1279  /* LP nonzeros are added as diagonal entries of the last block, the r-variable is variable nactivevars + 1 */
1280  sdpisolver->sdpa->inputElement((long long) (sdpisolver->nactivevars + 1), (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, 1.0, checkinput);
1281  }
1282  ++lhsrhscnt;
1283  }
1284 
1285  if ( rhspresent )
1286  {
1287  assert( 0 <= lastrow && lastrow < nlpcons );
1288 
1289  /* LP constraints are added as diagonal entries of the last block, right-hand side is added as variable zero with negative value */
1290  sdpisolver->sdpa->inputElement(0LL, (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, - lprhs[lastrow], checkinput);
1291 
1292  /* if we use a penalty formulation, add the r * Identity entry */
1293  if ( penaltyparam >= sdpisolver->epsilon )
1294  {
1295  /* LP nonzeros are added as diagonal entries of the last block (coming after the last SDP-block, with
1296  * blocks starting at 1, as are rows), the r-variable is variable nactivevars + 1 */
1297  sdpisolver->sdpa->inputElement((long long) (sdpisolver->nactivevars + 1), (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, 1.0, checkinput);
1298  }
1299  ++lhsrhscnt;
1300  }
1301 
1302  lastrow = lprow[i];
1303 
1304  /* determine presence of new lhs/rhs */
1305  if ( lplhs[lastrow] > - SCIPsdpiSolverInfinity(sdpisolver) )
1306  lhspresent = TRUE;
1307  else
1308  lhspresent = FALSE;
1309 
1310  if ( lprhs[lastrow] < SCIPsdpiSolverInfinity(sdpisolver) )
1311  rhspresent = TRUE;
1312  else
1313  rhspresent = FALSE;
1314  }
1315 
1316  /* if the variable is active */
1317  v = sdpisolver->inputtosdpamapper[lpcol[i]];
1318  if ( v > 0 )
1319  {
1320  /* add the LP-nonzero to the lhs-inequality if it exists: */
1321  if ( lhspresent )
1322  {
1323  /* LP nonzeros are added as diagonal entries of the last block */
1324  sdpisolver->sdpa->inputElement((long long) v, (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, lpval[i], checkinput);
1325  }
1326 
1327  /* add the LP-nonzero to the rhs-inequality if it exists: */
1328  if ( rhspresent )
1329  {
1330  /* LP nonzeros are added as diagonal entries of the last block; negative value because this is a <=-constraint, while SDPA works with >= */
1331  sdpisolver->sdpa->inputElement((long long) v, (long long) (nsdpablocks + 1), lhsrhscnt + lhspresent, lhsrhscnt + lhspresent, - lpval[i], checkinput);
1332  }
1333  }
1334  }
1335 
1336  /* take care of last row */
1337  if ( lhspresent )
1338  {
1339  assert( 0 <= lastrow && lastrow < nlpcons );
1340  sdpisolver->sdpa->inputElement(0LL, (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, lplhs[lastrow], checkinput);
1341  if ( penaltyparam >= sdpisolver->epsilon )
1342  sdpisolver->sdpa->inputElement((long long) (sdpisolver->nactivevars + 1), (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, 1.0, checkinput);
1343  ++lhsrhscnt;
1344  }
1345 
1346  if ( rhspresent )
1347  {
1348  assert( 0 <= lastrow && lastrow < nlpcons );
1349  sdpisolver->sdpa->inputElement(0LL, (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, - lprhs[lastrow], checkinput);
1350  if ( penaltyparam >= sdpisolver->epsilon )
1351  sdpisolver->sdpa->inputElement((long long) (sdpisolver->nactivevars + 1), (long long) (nsdpablocks + 1), lhsrhscnt, lhsrhscnt, 1.0, checkinput);
1352  ++lhsrhscnt;
1353  }
1354 
1355  assert( lhsrhscnt == (long long) nlpineqs + 1 );
1356  }
1357 
1358  /* inserting variable bounds */
1359  if ( sdpisolver->nvarbounds > 0 )
1360  {
1361  long long varboundcnt;
1362  int v;
1363 
1364  varboundcnt = (long long) nlpineqs + 1; /* +1 because of SDPA */
1365 
1366  for (i = 0; i < nvars; i++)
1367  {
1368  v = sdpisolver->inputtosdpamapper[i];
1369  if ( v > 0 )
1370  {
1371  assert( 0 < v && v <= sdpisolver->nactivevars );
1372  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, lb[i]) )
1373  {
1374  /* add bound as LP constraint for this variable */
1375  sdpisolver->sdpa->inputElement((long long) v, (long long) (nsdpablocks + 1), varboundcnt, varboundcnt, 1.0, checkinput);
1376 
1377  if ( REALABS(lb[i]) > sdpisolver->epsilon )
1378  {
1379  /* the bound is added as the rhs and therefore variable zero */
1380  sdpisolver->sdpa->inputElement(0LL, (long long) (nsdpablocks + 1), varboundcnt, varboundcnt, lb[i], checkinput);
1381  }
1382  ++varboundcnt;
1383  }
1384 
1385  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, ub[i]) )
1386  {
1387  /* add bound as LP constraint for this variable; - because of <= */
1388  sdpisolver->sdpa->inputElement((long long) v, (long long) (nsdpablocks + 1), varboundcnt, varboundcnt, -1.0, checkinput);
1389 
1390  if ( REALABS(ub[i]) > sdpisolver->epsilon )
1391  {
1392  /* the bound is added as the rhs and therefore variable zero, we multiply by -1 for <= */
1393  sdpisolver->sdpa->inputElement(0LL, (long long) (nsdpablocks + 1), varboundcnt, varboundcnt, -ub[i], checkinput);
1394  }
1395  ++varboundcnt;
1396  }
1397  }
1398  }
1399 
1400  if ( penaltyparam >= sdpisolver->epsilon && rbound )
1401  {
1402  /* we add the variable bound r >= 0 */
1403  sdpisolver->sdpa->inputElement((long long) (sdpisolver->nactivevars + 1), (long long) (nsdpablocks + 1), varboundcnt, varboundcnt, 1.0, checkinput);
1404  ++varboundcnt;
1405  }
1406  assert( varboundcnt == nlpineqs + 1 + sdpisolver->nvarbounds );
1407  }
1408 
1409 
1410  /* initialize settings (this needs to be done before inserting the problem as the initial point depends on the settings) */
1411  if ( penaltyparam >= sdpisolver->epsilon || startsettings == SCIP_SDPSOLVERSETTING_STABLE || startsettings == SCIP_SDPSOLVERSETTING_PENALTY )
1412  {
1413  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_STABLE_BUT_SLOW); /* if we already had problems with this problem, there is no reason to try fast */
1414  /* 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 */
1415  sdpisolver->sdpa->setParameterEpsilonStar(GAPTOLCHANGE * GAPTOLCHANGE * sdpisolver->gaptol);
1416  sdpisolver->sdpa->setParameterEpsilonDash(FEASTOLCHANGE * FEASTOLCHANGE * sdpisolver->sdpsolverfeastol);
1417  feastol = FEASTOLCHANGE * FEASTOLCHANGE * sdpisolver->sdpsolverfeastol;
1418  gaptol = sdpisolver->gaptol;
1419  SCIPdebugMessage("Start solving process with stable settings ...\n");
1420  }
1421  else if ( startsettings == SCIP_SDPSOLVERSETTING_UNSOLVED || startsettings == SCIP_SDPSOLVERSETTING_FAST )
1422  {
1423  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_UNSTABLE_BUT_FAST);
1424  sdpisolver->sdpa->setParameterEpsilonStar(sdpisolver->gaptol);
1425  sdpisolver->sdpa->setParameterEpsilonDash(sdpisolver->sdpsolverfeastol);
1426  feastol = sdpisolver->sdpsolverfeastol;
1427  gaptol = sdpisolver->gaptol;
1428  SCIPdebugMessage("Start solving process with fast settings\n");
1429  }
1430  else if ( startsettings == SCIP_SDPSOLVERSETTING_MEDIUM )
1431  {
1432  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_DEFAULT);
1433  /* as we want to solve with stable settings, we also update epsilon and the feasibility tolerance, as we skip the default settings, we multiply once */
1434  sdpisolver->sdpa->setParameterEpsilonStar(GAPTOLCHANGE * sdpisolver->gaptol);
1435  sdpisolver->sdpa->setParameterEpsilonDash(FEASTOLCHANGE * sdpisolver->sdpsolverfeastol);
1436  feastol = FEASTOLCHANGE * sdpisolver->sdpsolverfeastol;
1437  gaptol = sdpisolver->gaptol;
1438  SCIPdebugMessage("Start solving process with medium settings ...\n");
1439  }
1440  else
1441  {
1442  SCIPdebugMessage("Unknown setting for start-settings: %d!\n", startsettings); \
1443  return SCIP_LPERROR;
1444  }
1445 
1446  sdpisolver->sdpa->setParameterLowerBound(-1e20);
1447 
1448  /* set the objective limit */
1449  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
1450  sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
1451  else
1452  sdpisolver->sdpa->setParameterUpperBound(1e8);
1453 
1454  /* increase Lambda Star, this seems to help the numerics */
1455  sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
1456  //sdpisolver->sdpa->setParameterBetaStar(0.01);
1457  //sdpisolver->sdpa->setParameterBetaBar(0.02);
1458 
1459 #ifdef SCIP_MORE_DEBUG
1460  sdpisolver->sdpa->printParameters(stdout);
1461 #endif
1462 
1463  if ( sdpisolver->sdpinfo )
1464  sdpisolver->sdpa->setDisplay(stdout);
1465  else
1466  sdpisolver->sdpa->setDisplay(0);
1467 
1468 #ifdef SCIP_MORE_DEBUG
1469  FILE* fpOut = fopen("output.tmp", "w");
1470  if ( ! fpOut )
1471  exit(-1);
1472  sdpisolver->sdpa->setResultFile(fpOut);
1473 #endif
1474 
1475  /* if we want to use a starting point we have to tell SDPA to allocate memory for it */
1476  if ( starty != NULL && startZnblocknonz != NULL && startXnblocknonz != NULL && penaltyparam < sdpisolver->epsilon )
1477  sdpisolver->sdpa->setInitPoint(true);
1478  else
1479  sdpisolver->sdpa->setInitPoint(false);
1480 
1481  /* set the starting solution */
1482  if ( starty != NULL && startZnblocknonz != NULL && startXnblocknonz != NULL && penaltyparam < sdpisolver->epsilon )
1483  {
1484  int sdpaind;
1485  int varboundind;
1486 
1487  /* set dual vector */
1488  for (i = 1; i <= sdpisolver->nactivevars; i++)
1489  sdpisolver->sdpa->inputInitXVec((long long) i, starty[sdpisolver->sdpatoinputmapper[i - 1]]);/*lint !e747*/
1490 
1491  /* set SDP-blocks of dual matrix */
1492  for (b = 0; b < nsdpblocks; b++)
1493  {
1494  if ( blockindchanges[b] >= 0 )
1495  {
1496  for (i = 0; i < startZnblocknonz[b]; i++)
1497  {
1498  if ( indchanges[b][startZrow[b][i]] > -1 && indchanges[b][startZcol[b][i]] > -1 ) /* the row/col may have been fixed to zero in the meantime */
1499  {
1500  sdpisolver->sdpa->inputInitXMat((long long) b + 1 - blockindchanges[b], (long long) startZrow[b][i] + 1 - indchanges[b][startZrow[b][i]],
1501  (long long) startZcol[b][i] + 1 - indchanges[b][startZcol[b][i]], startZval[b][i]);
1502  }
1503  }
1504  }
1505  }
1506 
1507  /* set LP-block of dual matrix */
1508  for (i = 0; i < startZnblocknonz[nsdpblocks]; i++)
1509  {
1510  assert( startZrow[nsdpblocks][i] == startZcol[nsdpblocks][i] );
1511 
1512  if ( startZrow[nsdpblocks][i] < 2 * sdpisolver->nsdpalpcons ) /* linear constraint */
1513  {
1514  sdpaind = sdpisolver->rowmapper[startZrow[nsdpblocks][i]];
1515  if ( sdpaind > -1 )
1516  sdpisolver->sdpa->inputInitXMat((long long) nsdpblocks + 1, sdpaind, sdpaind, startZval[nsdpblocks][i]);
1517  }
1518  else /* varbound */
1519  {
1520  varboundind = sdpisolver->inputtovbmapper[startZrow[nsdpblocks][i] - 2 * sdpisolver->nsdpalpcons];
1521  if ( varboundind > -1 ) /* rhs */
1522  {
1523  sdpisolver->sdpa->inputInitXMat((long long) nsdpblocks + 1, nlpineqs + varboundind + 1,
1524  nlpineqs + varboundind + 1, startZval[nsdpblocks][i]);
1525  }
1526  }
1527  }
1528 
1529  /* set SDP-blocks of primal matrix */
1530  for (b = 0; b < nsdpblocks; b++)
1531  {
1532  if ( blockindchanges[b] >= 0 )
1533  {
1534  for (i = 0; i < startXnblocknonz[b]; i++)
1535  {
1536  if ( indchanges[b][startXrow[b][i]] > -1 && indchanges[b][startXcol[b][i]] > -1 ) /* the row/col may have been fixed to zero in the meantime */
1537  {
1538  sdpisolver->sdpa->inputInitYMat((long long) b + 1 - blockindchanges[b], (long long) startXrow[b][i] + 1 - indchanges[b][startXrow[b][i]],
1539  (long long) startXcol[b][i] + 1 - indchanges[b][startXcol[b][i]], startXval[b][i]);
1540  }
1541  }
1542  }
1543  }
1544 
1545  /* set LP-block of primal matrix */
1546  for (i = 0; i < startXnblocknonz[nsdpblocks]; i++)
1547  {
1548  assert( startXrow[nsdpblocks][i] == startXcol[nsdpblocks][i] );
1549 
1550  if ( startXrow[nsdpblocks][i] < 2 * sdpisolver->nsdpalpcons ) /* linear constraint */
1551  {
1552  sdpaind = sdpisolver->rowmapper[startXrow[nsdpblocks][i]];
1553  if ( sdpaind > -1 )
1554  sdpisolver->sdpa->inputInitYMat((long long) nsdpblocks + 1, sdpaind, sdpaind, startXval[nsdpblocks][i]);
1555  }
1556  else /* varbound */
1557  {
1558  varboundind = sdpisolver->inputtovbmapper[startXrow[nsdpblocks][i] - 2 * sdpisolver->nsdpalpcons];
1559  if ( varboundind > -1 ) /*lhs */
1560  {
1561  sdpisolver->sdpa->inputInitYMat((long long) nsdpblocks + 1, nlpineqs + varboundind + 1,
1562  nlpineqs + varboundind + 1, startXval[nsdpblocks][i]);
1563  }
1564  }
1565  }
1566  }
1567  else if ( penaltyparam >= sdpisolver->epsilon )
1568  SCIPdebugMessage("Skipping insertion of starting point for penalty formulation.\n");
1569 
1570  /* transform the matrices to a more efficient form */
1571  sdpisolver->sdpa->initializeUpperTriangle(checkinput);
1572  sdpisolver->sdpa->initializeSolve();
1573 
1574 #ifdef SCIP_DEBUG_PRINTTOFILE
1575  /* if necessary, dump input data and initial point */
1576  sdpisolver->sdpa->writeInputSparse(const_cast<char*>("sdpa.dat-s"), const_cast<char*>("%+8.3e"));
1577  sdpisolver->sdpa->writeInitSparse(const_cast<char*>("sdpa.ini-s"), const_cast<char*>("%+8.3e"));
1578 #endif
1579 
1580  SCIPdebugMessage("Calling SDPA solve (SDP: %d) ...\n", sdpisolver->sdpcounter);
1581 
1582  /* check if we want to save a preoptimal solution for warmstarting purposes */
1583  if ( usepreoptimalsol && (startsettings == SCIP_SDPSOLVERSETTING_UNSOLVED || startsettings == SCIP_SDPSOLVERSETTING_FAST) )
1584  {
1585  SCIP_Real* sdpasol;
1586  int sdpablocksize;
1587 
1588  /* first solve up to gaptol */
1589  sdpisolver->sdpa->setParameterEpsilonStar(sdpisolver->preoptimalgap);
1590  sdpisolver->sdpa->setParameterEpsilonDash(sdpisolver->preoptimalgap);
1591  sdpisolver->sdpa->solve();
1592 
1593  /* save preoptimal solution (only if solving succeeded) */
1594  if ( SCIPsdpiSolverIsAcceptable(sdpisolver) )
1595  {
1596  SCIPdebugMessage("Saving preoptimal solution for warmstarting purposes\n");
1597  sdpasol = sdpisolver->sdpa->getResultXVec();
1598 
1599  /* copy dual solution vector */
1600  for (i = 0; i < sdpisolver->nactivevars; i++)
1601  sdpisolver->preoptimalsol[i] = sdpasol[i];
1602 
1603  /* copy primal matrix */
1604  for (b = 0; b < sdpisolver->sdpa->getBlockNumber(); b++)
1605  {
1606  sdpasol = sdpisolver->sdpa->getResultYMat(b + 1);
1607  sdpablocksize = (int) sdpisolver->sdpa->getBlockSize(b + 1);
1608  assert( sdpablocksize <= 2 * (sdpisolver->maxnlpcons + sdpisolver->maxnvars) );
1609  if ( sdpisolver->sdpa->getBlockType(b + 1) == SDPA::LP )
1610  {
1611  /* only the diagonal entries need to be saved */
1612  for (i = 0; i < sdpablocksize; i++)
1613  sdpisolver->preoptimalsolxlp[i] = sdpasol[i];
1614  }
1615  else
1616  {
1617  /* TODO: think about saving this in sparse format or at least only the upper triangular matrix */
1618  for (i = 0; i < sdpablocksize * sdpablocksize; i++)
1619  sdpisolver->preoptimalsolx[b][i] = sdpasol[i];
1620  }
1621  }
1622 
1623  sdpisolver->preoptimalsolexists = TRUE;
1624  }
1625  else
1626  {
1627  SCIPdebugMessage("Solving to gaptol failed! No preoptimal solution available.\n");
1628  }
1629 
1630  /* add number of iterations */
1631  sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1632 
1633  /* copy current iterate and resolve with real tolerance */
1634  sdpisolver->sdpa->setInitPoint(true);
1635  sdpisolver->sdpa->copyCurrentToInit();
1636  sdpisolver->sdpa->setParameterEpsilonStar(sdpisolver->gaptol);
1637  sdpisolver->sdpa->setParameterEpsilonDash(sdpisolver->sdpsolverfeastol);
1638  sdpisolver->sdpa->solve();
1639  }
1640  else
1641  sdpisolver->sdpa->solve();
1642 
1643  sdpisolver->solved = TRUE;
1644 
1645  /* update number of SDP-iterations and -calls */
1646  sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1647  ++sdpisolver->nsdpcalls;
1648 
1649 #ifdef SCIP_DEBUG
1650  /* print the phase value , i.e. whether solving was successfull */
1651  sdpisolver->sdpa->getPhaseString((char*)phase_string);
1652  SCIPdebugMessage("SDPA solving finished with status %s (primal and dual here are switched in contrast to our formulation)\n", phase_string);
1653 #endif
1654 
1655  /* remember settings */
1656  if ( SCIPsdpiSolverIsAcceptable(sdpisolver) && penaltyparam < sdpisolver->epsilon )
1657  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_FAST;
1658  else if ( penaltyparam >= sdpisolver->epsilon )
1659  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_PENALTY;
1660 
1661  /* if the problem has been stably solved but did not reach the required feasibility tolerance, even though the solver
1662  * reports feasibility, resolve it with adjusted tolerance */
1663  SCIP_CALL( checkFeastolAndResolve(sdpisolver, penaltyparam, nvars, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
1664  sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
1665  nremovedinds, blockindchanges, nlpcons, lplhs, lprhs, lpnnonz, lprow, lpcol, lpval, &feastol, &gaptol) );
1666 
1667  /* 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
1668  * again with more stable parameters */
1669  if ( ! SCIPsdpiSolverIsAcceptable(sdpisolver) && penaltyparam < sdpisolver->epsilon &&
1670  (startsettings == SCIP_SDPSOLVERSETTING_UNSOLVED || startsettings == SCIP_SDPSOLVERSETTING_FAST) )
1671  {
1672  SCIPdebugMessage("Numerical troubles -- solving SDP %d again ...\n", sdpisolver->sdpcounter);
1673 
1674  /* initialize settings */
1675  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_DEFAULT);
1676  sdpisolver->sdpa->setParameterEpsilonStar(GAPTOLCHANGE * sdpisolver->gaptol);
1677  sdpisolver->sdpa->setParameterEpsilonDash(FEASTOLCHANGE * sdpisolver->sdpsolverfeastol);
1678  sdpisolver->sdpa->setParameterLowerBound(-1e20);
1679 
1680  /* set the objective limit */
1681  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
1682  sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
1683  else
1684  sdpisolver->sdpa->setParameterUpperBound(1e8);
1685 
1686  /* increase Lambda Star, this seems to help the numerics */
1687  sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
1688 
1689 #ifdef SCIP_MORE_DEBUG
1690  sdpisolver->sdpa->printParameters(stdout);
1691 #endif
1692  sdpisolver->sdpa->setInitPoint(false);
1693 #ifdef SDPA_RESETPARAMS
1694  sdpisolver->sdpa->resetParameters();
1695 #else
1696  sdpisolver->sdpa->initializeSolve();
1697 #endif
1698  sdpisolver->sdpa->solve();
1699  sdpisolver->solved = TRUE;
1700 
1701  /* update number of SDP-iterations and -calls */
1702  sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1703  ++sdpisolver->nsdpcalls;
1704 
1705  /* remember setting */
1706  if ( SCIPsdpiSolverIsAcceptable(sdpisolver) )
1707  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_MEDIUM;
1708 
1709 #ifdef SCIP_DEBUG
1710  /* print the phase value , i.e. whether solving was successfull */
1711  sdpisolver->sdpa->getPhaseString((char*)phase_string);
1712  SCIPdebugMessage("SDPA solving finished with status %s (primal and dual here are switched in contrast to our formulation)\n", phase_string);
1713 #endif
1714 
1715  /* if the problem has been stably solved but did not reach the required feasibility tolerance, even though the solver
1716  * reports feasibility, resolve it with adjusted tolerance */
1717  SCIP_CALL( checkFeastolAndResolve(sdpisolver, penaltyparam, nvars, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
1718  sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
1719  nremovedinds, blockindchanges, nlpcons, lplhs, lprhs, lpnnonz, lprow, lpcol, lpval, &feastol, &gaptol) );
1720  }
1721 
1722  /* if we still didn't converge, and did not yet use the stable settings, set the parameters even more conservativly */
1723  if ( (! SCIPsdpiSolverIsAcceptable(sdpisolver)) && penaltyparam < sdpisolver->epsilon &&
1724  (startsettings == SCIP_SDPSOLVERSETTING_UNSOLVED || startsettings == SCIP_SDPSOLVERSETTING_FAST || startsettings == SCIP_SDPSOLVERSETTING_MEDIUM) )
1725  {
1726  SCIPdebugMessage("Numerical troubles -- solving SDP %d again^2 ...\n", sdpisolver->sdpcounter);
1727 
1728  /* initialize settings */
1729  sdpisolver->sdpa->setParameterType(SDPA::PARAMETER_STABLE_BUT_SLOW);
1730  sdpisolver->sdpa->setParameterEpsilonStar(GAPTOLCHANGE * GAPTOLCHANGE * sdpisolver->gaptol);
1731  sdpisolver->sdpa->setParameterEpsilonDash(FEASTOLCHANGE * FEASTOLCHANGE * sdpisolver->sdpsolverfeastol);
1732  sdpisolver->sdpa->setParameterLowerBound(-1e20);
1733 
1734  /* set the objective limit */
1735  if ( ! SCIPsdpiSolverIsInfinity(sdpisolver, sdpisolver->objlimit) )
1736  sdpisolver->sdpa->setParameterUpperBound(sdpisolver->objlimit);
1737  else
1738  sdpisolver->sdpa->setParameterUpperBound(1e8);
1739 
1740  /* increase Lambda Star, this seems to help the numerics */
1741  sdpisolver->sdpa->setParameterLambdaStar(sdpisolver->lambdastar);
1742 
1743 #ifdef SCIP_MORE_DEBUG
1744  sdpisolver->sdpa->printParameters(stdout);
1745 #endif
1746  sdpisolver->sdpa->setInitPoint(false);
1747 #ifdef SDPA_RESETPARAMS
1748  sdpisolver->sdpa->resetParameters();
1749 #else
1750  sdpisolver->sdpa->initializeSolve();
1751 #endif
1752  sdpisolver->sdpa->solve();
1753  sdpisolver->solved = TRUE;
1754 
1755  /* update number of SDP-iterations and -calls */
1756  sdpisolver->niterations += (int) sdpisolver->sdpa->getIteration();
1757  ++sdpisolver->nsdpcalls;
1758 
1759  /* remember setting */
1760  if ( SCIPsdpiSolverIsAcceptable(sdpisolver) )
1761  sdpisolver->usedsetting = SCIP_SDPSOLVERSETTING_STABLE;
1762 
1763 #ifdef SCIP_DEBUG
1764  /* print the phase value , i.e. whether solving was successfull */
1765  sdpisolver->sdpa->getPhaseString((char*)phase_string);
1766  SCIPdebugMessage("SDPA solving finished with status %s (primal and dual here are switched in constrast to our formulation)\n", phase_string);
1767 #endif
1768 
1769  /* if the problem has been stably solved but did not reach the required feasibility tolerance, even though the solver
1770  * reports feasibility, resolve it with adjusted tolerance */
1771  SCIP_CALL( checkFeastolAndResolve(sdpisolver, penaltyparam, nvars, lb, ub, nsdpblocks, sdpblocksizes, sdpnblockvars, sdpconstnnonz,
1772  sdpconstnblocknonz, sdpconstrow, sdpconstcol, sdpconstval, sdpnnonz, sdpnblockvarnonz, sdpvar, sdprow, sdpcol, sdpval, indchanges,
1773  nremovedinds, blockindchanges, nlpcons, lplhs, lprhs, lpnnonz, lprow, lpcol, lpval, &feastol, &gaptol) );
1774  }
1775 
1776 #ifdef SCIP_MORE_DEBUG
1777  (void) fclose(fpOut);
1778 #endif
1779 
1780  /* if we solved a penalty formulation, check if the solution is feasible for the original problem (which is the case iff r < feastol) */
1781  if ( penaltyparam >= sdpisolver->epsilon )
1782  {
1783  SCIP_Real* sdpasol;
1784  SCIP_Real* X;
1785  int nblockssdpa;
1786  int nrow;
1787  SCIP_Real trace = 0.0;
1788 
1789  /* in the second case we have r as an additional variable */
1790  assert( (sdpisolver->nactivevars + 1 == sdpisolver->sdpa->getConstraintNumber()) ); /*lint !e776*/
1791 
1792  sdpasol = sdpisolver->sdpa->getResultXVec();
1793 
1794  /* we get r as the last variable in SDPA */
1795  assert( feasorig != NULL );
1796  if ( sdpasol[sdpisolver->nactivevars] < sdpisolver->feastol )
1797  *feasorig = TRUE;
1798  else
1799  *feasorig = FALSE;
1800 
1801  /* only set sdpisolver->feasorig to true if we solved with objective, because only in this case we want to compute
1802  * the objective value by hand since it is numerically more stable then the result returned by SDPA */
1803  if ( withobj )
1804  sdpisolver->feasorig = *feasorig;
1805 
1806  /* if r > 0 or we are in debug mode, also check the primal bound */
1807 #ifdef NDEBUG
1808  if ( ! *feasorig && penaltybound != NULL )
1809  {
1810 #endif
1811  SCIPdebugMessage("Solution not feasible in original problem, r = %f\n", sdpasol[sdpisolver->nactivevars]);
1812 
1813  /* compute Tr(X) */
1814 
1815  /* iterate over all blocks (SDPA starts counting at one and includes the LP block) */
1816  nblockssdpa = (int) sdpisolver->sdpa->getBlockNumber();
1817  for (b = 1; b <= nblockssdpa; b++)
1818  {
1819  /* get the block from SDPA */
1820  X = sdpisolver->sdpa->getResultYMat((long long) b);/*lint !e747*/
1821  nrow = (int) sdpisolver->sdpa->getBlockSize((long long) b);/*lint !e747*/
1822  assert( nrow >= 0 );
1823 
1824  /* if it is the LP-block, we omit the variable bounds as the penalty variable is not added to them */
1825  if ( sdpisolver->sdpa->getBlockType((long long) b) == SDPA::LP )/*lint !e747*/
1826  {
1827  /* iterate over all diagonal entries (until we reach the varbound part), adding them to the trace */
1828  for (i = 0; i < nrow - sdpisolver->nvarbounds; i++)
1829  trace += X[i]; /* get entry (i+1,i+1) for the diagonal matrix X */
1830  }
1831  else
1832  {
1833  /* iterate over all diagonal entries and add them to the trace */
1834  for (i = 0; i < nrow; i++)
1835  trace += X[i + i*nrow];/*lint !e679*/ /* get entry (i+1,i+1) in X */
1836  }
1837  }
1838 
1839  /* if the relative gap is smaller than the tolerance, we return equality */
1840  if ( (penaltyparam - trace) / penaltyparam < PENALTYBOUNDTOL )/*lint !e414*/
1841  {
1842  if ( penaltybound != NULL )
1843  *penaltybound = TRUE;
1844 
1845  SCIPdebugMessage("Tr(X) = %f == %f = Gamma, penalty formulation not exact, Gamma should be increased or problem is infeasible\n",
1846  trace, penaltyparam);
1847  }
1848  else if ( penaltybound != NULL )
1849  *penaltybound = FALSE;
1850 
1851 #ifdef NDEBUG
1852  }
1853 #endif
1854  }
1855  return SCIP_OKAY;
1856 }
1857 
1863 /*
1864  * Solution Information Methods
1865  */
1866 
1871 SCIP_Bool SCIPsdpiSolverWasSolved(
1872  SCIP_SDPISOLVER* sdpisolver
1873  )
1874 {/*lint !e1784*/
1875  assert( sdpisolver != NULL );
1876  return sdpisolver->solved;
1877 }
1878 
1885 SCIP_Bool SCIPsdpiSolverFeasibilityKnown(
1886  SCIP_SDPISOLVER* sdpisolver
1887  )
1888 {/*lint !e1784*/
1889  SDPA::PhaseType phasetype;
1890 
1891  assert( sdpisolver != NULL );
1892  assert( sdpisolver->sdpa != NULL);
1893  CHECK_IF_SOLVED_BOOL( sdpisolver );
1894 
1895  phasetype = sdpisolver->sdpa->getPhaseValue();
1896 
1897  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
1898  return FALSE;
1899 
1900  return TRUE;
1901 }
1902 
1904 SCIP_RETCODE SCIPsdpiSolverGetSolFeasibility(
1905  SCIP_SDPISOLVER* sdpisolver,
1906  SCIP_Bool* primalfeasible,
1907  SCIP_Bool* dualfeasible
1908  )
1909 {/*lint !e1784*/
1910  SDPA::PhaseType phasetype;
1911 
1912  assert( sdpisolver != NULL );
1913  assert( sdpisolver->sdpa != NULL );
1914  assert( primalfeasible != NULL );
1915  assert( dualfeasible != NULL );
1916  CHECK_IF_SOLVED( sdpisolver );
1917 
1918  phasetype = sdpisolver->sdpa->getPhaseValue();
1919 
1920  switch ( phasetype )/*lint --e{788}*/
1921  {
1922  case SDPA::pdOPT:
1923  *primalfeasible = TRUE;
1924  *dualfeasible = TRUE;
1925  break;
1926  case SDPA::pdFEAS:
1927  *primalfeasible = TRUE;
1928  *dualfeasible = TRUE;
1929  break;
1930  case SDPA::pFEAS_dINF:
1931  *primalfeasible = TRUE;
1932  *dualfeasible = FALSE;
1933  break;
1934  case SDPA::pINF_dFEAS:
1935  *primalfeasible = FALSE;
1936  *dualfeasible = TRUE;
1937  break;
1938  case SDPA::pUNBD:
1939  *primalfeasible = TRUE;
1940  *dualfeasible = FALSE;
1941  SCIPdebugMessage("SDPA stopped because dual objective became smaller than lower bound\n");
1942  break;
1943  case SDPA::dUNBD:
1944  *primalfeasible = FALSE;
1945  *dualfeasible = TRUE;
1946  SCIPdebugMessage("SDPA stopped because primal objective became bigger than upper bound\n");
1947  break;
1948  default: /* contains noInfo, pFeas, dFeas, pdInf */
1949  SCIPerrorMessage("SDPA doesn't know if primal and dual solutions are feasible\n");
1950  return SCIP_LPERROR;
1951  }
1952 
1953  return SCIP_OKAY;
1954 }
1955 
1959 SCIP_Bool SCIPsdpiSolverIsPrimalUnbounded(
1960  SCIP_SDPISOLVER* sdpisolver
1961  )
1962 {/*lint !e1784*/
1963  SDPA::PhaseType phasetype;
1964 
1965  assert( sdpisolver != NULL );
1966  assert( sdpisolver->sdpa != NULL );
1967  CHECK_IF_SOLVED_BOOL( sdpisolver );
1968 
1969  phasetype = sdpisolver->sdpa->getPhaseValue();
1970 
1971  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::pdINF )
1972  {
1973  SCIPdebugMessage("SDPA doesn't know if primal problem is unbounded");
1974  return FALSE;
1975  }
1976  else if ( phasetype == SDPA::pFEAS_dINF )
1977  return TRUE;
1978  else if ( phasetype == SDPA::pUNBD )
1979  {
1980  SCIPdebugMessage("SDPA was stopped because primal objective became bigger than upper bound");
1981  return TRUE;
1982  }
1983 
1984  return FALSE;
1985 }
1986 
1990 SCIP_Bool SCIPsdpiSolverIsPrimalInfeasible(
1991  SCIP_SDPISOLVER* sdpisolver
1992  )
1993 {/*lint !e1784*/
1994  SDPA::PhaseType phasetype;
1995 
1996  assert( sdpisolver != NULL );
1997  assert( sdpisolver->sdpa != NULL );
1998  CHECK_IF_SOLVED_BOOL( sdpisolver );
1999 
2000  phasetype = sdpisolver->sdpa->getPhaseValue();
2001 
2002  if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
2003  {
2004  SCIPdebugMessage("SDPA doesn't know if primal problem is infeasible");
2005  return FALSE;
2006  }
2007  else if ( phasetype == SDPA::pINF_dFEAS )
2008  return TRUE;
2009  else if ( phasetype == SDPA::dUNBD )
2010  {
2011  SCIPdebugMessage("SDPA was stopped because dual objective became smaller than lower bound");
2012  return TRUE;
2013  }
2014 
2015  return FALSE;
2016 }
2017 
2021 SCIP_Bool SCIPsdpiSolverIsPrimalFeasible(
2022  SCIP_SDPISOLVER* sdpisolver
2023  )
2024 {/*lint !e1784*/
2025  SDPA::PhaseType phasetype;
2026 
2027  assert( sdpisolver != NULL );
2028  assert( sdpisolver->sdpa != NULL );
2029  CHECK_IF_SOLVED_BOOL( sdpisolver );
2030 
2031  phasetype = sdpisolver->sdpa->getPhaseValue();
2032 
2033  if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
2034  {
2035  SCIPdebugMessage("SDPA doesn't know if primal problem is feasible");
2036  return FALSE;
2037  }
2038  else if ( phasetype == SDPA::pFEAS_dINF || phasetype == SDPA::pdOPT || phasetype == SDPA::pFEAS || phasetype == SDPA::pdFEAS )
2039  return TRUE;
2040  else if ( phasetype == SDPA::pUNBD )
2041  {
2042  SCIPdebugMessage("SDPA was stopped because dual objective became smaller than lower bound");
2043  return TRUE;
2044  }
2045 
2046  return FALSE;
2047 }
2048 
2052 SCIP_Bool SCIPsdpiSolverIsDualUnbounded(
2053  SCIP_SDPISOLVER* sdpisolver
2054  )
2055 {/*lint !e1784*/
2056  SDPA::PhaseType phasetype;
2057 
2058  assert( sdpisolver != NULL );
2059  assert( sdpisolver->sdpa != NULL);
2060  CHECK_IF_SOLVED_BOOL( sdpisolver );
2061 
2062  phasetype = sdpisolver->sdpa->getPhaseValue();
2063 
2064  if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
2065  {
2066  SCIPdebugMessage("SDPA doesn't know if dual problem is unbounded");
2067  return FALSE;
2068  }
2069  else if ( phasetype == SDPA::pINF_dFEAS )
2070  return TRUE;
2071  else if ( phasetype == SDPA::dUNBD )
2072  {
2073  SCIPdebugMessage("SDPA was stopped because dual objective became smaller than lower bound");
2074  return TRUE;
2075  }
2076 
2077  return FALSE;
2078 }
2079 
2083 SCIP_Bool SCIPsdpiSolverIsDualInfeasible(
2084  SCIP_SDPISOLVER* sdpisolver
2085  )
2086 {/*lint !e1784*/
2087  SDPA::PhaseType phasetype;
2088 
2089  assert( sdpisolver != NULL );
2090  assert( sdpisolver->sdpa != NULL);
2091  CHECK_IF_SOLVED_BOOL( sdpisolver );
2092 
2093  phasetype = sdpisolver->sdpa->getPhaseValue();
2094 
2095  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::pdINF )
2096  {
2097  SCIPdebugMessage("SDPA doesn't know if dual problem is infeasible.\n");
2098  return FALSE;
2099  }
2100  else if ( phasetype == SDPA::pFEAS_dINF )
2101  return TRUE;
2102  else if ( phasetype == SDPA::pUNBD )
2103  {
2104  SCIPdebugMessage("SDPA was stopped because primal objective became bigger than upper bound");
2105  return TRUE;
2106  }
2107 
2108  return FALSE;
2109 }
2110 
2114 SCIP_Bool SCIPsdpiSolverIsDualFeasible(
2115  SCIP_SDPISOLVER* sdpisolver
2116  )
2117 {/*lint !e1784*/
2118  SDPA::PhaseType phasetype;
2119 
2120  assert( sdpisolver != NULL );
2121  assert( sdpisolver->sdpa != NULL);
2122  CHECK_IF_SOLVED_BOOL( sdpisolver );
2123 
2124  phasetype = sdpisolver->sdpa->getPhaseValue();
2125 
2126  if ( phasetype == SDPA::noINFO || phasetype == SDPA::dFEAS || phasetype == SDPA::pdINF )
2127  {
2128  SCIPdebugMessage("SDPA doesn't know if primal problem is feasible");
2129  return FALSE;
2130  }
2131  else if ( phasetype == SDPA::pINF_dFEAS || phasetype == SDPA::pdOPT || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
2132  return TRUE;
2133  else if ( phasetype == SDPA::dUNBD )
2134  {
2135  SCIPdebugMessage("SDPA was stopped because dual objective became smaller than lower bound");
2136  return TRUE;
2137  }
2138 
2139  return FALSE;
2140 }
2141 
2143 SCIP_Bool SCIPsdpiSolverIsConverged(
2144  SCIP_SDPISOLVER* sdpisolver
2145  )
2146 {/*lint !e1784*/
2147  SDPA::PhaseType phasetype;
2148 
2149  assert( sdpisolver != NULL );
2150  assert( sdpisolver->sdpa != NULL );
2151  CHECK_IF_SOLVED_BOOL( sdpisolver );
2152 
2153  phasetype = sdpisolver->sdpa->getPhaseValue();
2154 
2155  if ( phasetype == SDPA::pdOPT )
2156  return TRUE;
2157 
2158  return FALSE;
2159 }
2160 
2162 SCIP_Bool SCIPsdpiSolverIsObjlimExc(
2163  SCIP_SDPISOLVER* sdpisolver
2164  )
2165 {/*lint !e1784*/
2166  SDPA::PhaseType phasetype;
2167 
2168  assert( sdpisolver != NULL );
2169  assert( sdpisolver->sdpa != NULL );
2170  CHECK_IF_SOLVED_BOOL( sdpisolver );
2171 
2172  phasetype = sdpisolver->sdpa->getPhaseValue();
2173 
2174  if ( phasetype == SDPA::pUNBD )
2175  return TRUE;
2176 
2177  return FALSE;
2178 }
2179 
2181 SCIP_Bool SCIPsdpiSolverIsIterlimExc(
2182  SCIP_SDPISOLVER* sdpisolver
2183  )
2184 {/*lint !e1784*/
2185  SDPA::PhaseType phasetype;
2186 
2187  assert( sdpisolver != NULL );
2188  assert( sdpisolver->sdpa != NULL);
2189  CHECK_IF_SOLVED_BOOL( sdpisolver );
2190 
2191  phasetype = sdpisolver->sdpa->getPhaseValue();
2192 
2193  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
2194  {
2195  if ( sdpisolver->sdpa->getParameterMaxIteration() == sdpisolver->sdpa->getIteration() )
2196  return TRUE;
2197  }
2198 
2199  return FALSE;
2200 }
2201 
2203 SCIP_Bool SCIPsdpiSolverIsTimelimExc(
2204  SCIP_SDPISOLVER* sdpisolver
2205  )
2206 {/*lint !e1784*/
2207  assert( sdpisolver != NULL );
2208 
2209  return sdpisolver->timelimit;
2210 }
2211 
2223 int SCIPsdpiSolverGetInternalStatus(
2224  SCIP_SDPISOLVER* sdpisolver
2225  )
2226 {/*lint !e1784*/
2227  SDPA::PhaseType phasetype;
2228 
2229  assert( sdpisolver != NULL );
2230  assert( sdpisolver->sdpa != NULL );
2231 
2232  if ( sdpisolver->sdpa == NULL || (! sdpisolver->solved) )
2233  return -1;
2234 
2235  phasetype = sdpisolver->sdpa->getPhaseValue();
2236 
2237  if ( phasetype == SDPA::pdOPT || phasetype == SDPA::pFEAS_dINF || phasetype == SDPA::pINF_dFEAS )
2238  return 0;
2239  if ( phasetype == SDPA::pdINF )
2240  return 1;
2241  if ( phasetype == SDPA::pUNBD )
2242  return 3;
2243  if ( phasetype == SDPA::noINFO || phasetype == SDPA::pFEAS || phasetype == SDPA::dFEAS || phasetype == SDPA::pdFEAS )
2244  return 4;
2245  else /* should include dUNBD */
2246  return 7;
2247 }
2248 
2250 SCIP_Bool SCIPsdpiSolverIsOptimal(
2251  SCIP_SDPISOLVER* sdpisolver
2252  )
2253 {/*lint !e1784*/
2254  SDPA::PhaseType phasetype;
2255 
2256  assert( sdpisolver != NULL );
2257  assert( sdpisolver->sdpa != NULL );
2258 
2259  if ( ! sdpisolver->solved )
2260  return FALSE;
2261 
2262  phasetype = sdpisolver->sdpa->getPhaseValue();
2263 
2264  if ( phasetype == SDPA::pdOPT )
2265  return TRUE;
2266 
2267  return FALSE;
2268 }
2269 
2273 SCIP_Bool SCIPsdpiSolverIsAcceptable(
2274  SCIP_SDPISOLVER* sdpisolver
2275  )
2276 {/*lint !e1784*/
2277  SDPA::PhaseType phasetype;
2278 
2279  assert( sdpisolver != NULL );
2280  assert( sdpisolver->sdpa != NULL );
2281 
2282  if ( sdpisolver->timelimit )
2283  return FALSE;
2284 
2285  if ( ! sdpisolver->solved )
2286  return FALSE;
2287 
2288  phasetype = sdpisolver->sdpa->getPhaseValue();
2289 
2290  /* we are happy if we converged, or we reached the objective limit (pUNBD) or we could show that our problem is
2291  * infeasible [except for numerics], or unbounded */
2292  if ( SCIPsdpiSolverIsConverged(sdpisolver) || phasetype == SDPA::pUNBD || phasetype == SDPA::pINF_dFEAS || phasetype == SDPA::pFEAS_dINF )
2293  return TRUE;
2294 
2295  return FALSE;
2296 }
2297 
2299 SCIP_RETCODE SCIPsdpiSolverIgnoreInstability(
2300  SCIP_SDPISOLVER* sdpisolver,
2301  SCIP_Bool* success
2302  )
2303 {/*lint !e1784*/
2304  SCIPdebugMessage("Not implemented yet\n");
2305 
2306  return SCIP_LPERROR;
2307 }/*lint !e715*/
2308 
2310 SCIP_RETCODE SCIPsdpiSolverGetObjval(
2311  SCIP_SDPISOLVER* sdpisolver,
2312  SCIP_Real* objval
2313  )
2314 {/*lint !e1784*/
2315  assert( sdpisolver != NULL );
2316  assert( sdpisolver->sdpa != NULL );
2317  assert( objval != NULL );
2318  CHECK_IF_SOLVED( sdpisolver );
2319 
2320  if ( sdpisolver->penalty && ( ! sdpisolver->feasorig ) )
2321  {
2322  *objval = sdpisolver->sdpa->getPrimalObj();
2323 #ifndef NDEBUG
2324  SCIP_Real primalval = sdpisolver->sdpa->getDualObj();
2325  SCIP_Real gap = (REALABS(*objval - primalval) / (0.5 * (REALABS(primalval) + REALABS(*objval)))); /* duality gap used in SDPA */
2326  if ( gap > sdpisolver->gaptol )
2327  {
2328  SCIPdebugMessage("Attention: got objective value (before adding values of fixed variables) of %f in SCIPsdpiSolverGetSol, "
2329  "but primal objective is %g with duality gap %g!\n", *objval, primalval, gap );
2330  }
2331 #endif
2332  }
2333  else
2334  {
2335  SCIP_Real* sdpasol;
2336  int v;
2337 
2338  /* since the objective value given by SDPA sometimes differs slightly from the correct value for the given solution,
2339  * we get the solution from SDPA and compute the correct objective value */
2340  assert( (sdpisolver->penalty && sdpisolver->nactivevars + 1 == sdpisolver->sdpa->getConstraintNumber()) || /*lint !e776*/
2341  sdpisolver->nactivevars == sdpisolver->sdpa->getConstraintNumber() ); /* in the second case we have r as an additional variable */
2342  sdpasol = sdpisolver->sdpa->getResultXVec();
2343 
2344  *objval = 0.0;
2345  for (v = 0; v < sdpisolver->nactivevars; v++)
2346  *objval += sdpasol[v] * sdpisolver->objcoefs[v];
2347  }
2348 
2349  /* as we didn't add the fixed (lb = ub) variables to sdpa, we have to add their contributions to the objective by hand */
2350  *objval += sdpisolver->fixedvarsobjcontr;
2351 
2352  return SCIP_OKAY;
2353 }
2354 
2359 SCIP_RETCODE SCIPsdpiSolverGetSol(
2360  SCIP_SDPISOLVER* sdpisolver,
2361  SCIP_Real* objval,
2362  SCIP_Real* dualsol,
2363  int* dualsollength
2365  )
2366 {/*lint !e1784*/
2367  SCIP_Real* sdpasol;
2368  int v;
2369 
2370  assert( sdpisolver != NULL );
2371  assert( sdpisolver->sdpa != NULL );
2372  assert( dualsollength != NULL );
2373  CHECK_IF_SOLVED( sdpisolver );
2374 
2375  if ( objval != NULL )
2376  {
2377  if ( sdpisolver->penalty && ! sdpisolver->feasorig )
2378  {
2379  *objval = sdpisolver->sdpa->getPrimalObj();
2380 #ifndef NDEBUG
2381  SCIP_Real primalval = sdpisolver->sdpa->getDualObj();
2382  SCIP_Real gap = (REALABS(*objval - primalval) / (0.5 * (REALABS(primalval) + REALABS(*objval)))); /* duality gap used in SDPA */
2383  if ( gap > sdpisolver->gaptol )
2384  {
2385  SCIPdebugMessage("Attention: got objective value (before adding values of fixed variables) of %f in SCIPsdpiSolverGetSol, "
2386  "but primal objective is %f with duality gap %f!\n", *objval, primalval, gap );
2387  }
2388 #endif
2389  }
2390  else
2391  {
2392  /* since the objective value given by SDPA sometimes differs slightly from the correct value for the given solution,
2393  * we get the solution from SDPA and compute the correct objective value */
2394  assert( (sdpisolver->penalty && sdpisolver->nactivevars + 1 == sdpisolver->sdpa->getConstraintNumber()) || /*lint !e776*/
2395  sdpisolver->nactivevars == sdpisolver->sdpa->getConstraintNumber() ); /* in the second case we have r as an additional variable */
2396  sdpasol = sdpisolver->sdpa->getResultXVec();
2397 
2398  *objval = 0.0;
2399  for (v = 0; v < sdpisolver->nactivevars; v++)
2400  *objval += sdpasol[v] * sdpisolver->objcoefs[v];
2401  }
2402 
2403  /* as we didn't add the fixed (lb = ub) variables to sdpa, we have to add their contributions to the objective by hand */
2404  *objval += sdpisolver->fixedvarsobjcontr;
2405  }
2406 
2407  if ( *dualsollength > 0 )
2408  {
2409  assert( dualsol != NULL );
2410  if ( *dualsollength < sdpisolver->nvars )
2411  {
2412  SCIPdebugMessage("The given array in SCIPsdpiSolverGetSol only had length %d, but %d was needed", *dualsollength, sdpisolver->nvars);
2413  *dualsollength = sdpisolver->nvars;
2414 
2415  return SCIP_OKAY;
2416  }
2417 
2418  /* get the solution from sdpa */
2419  assert( (sdpisolver->penalty && sdpisolver->nactivevars + 1 == sdpisolver->sdpa->getConstraintNumber()) || /*lint !e776*/
2420  sdpisolver->nactivevars == sdpisolver->sdpa->getConstraintNumber() ); /* in the second case we have r as an additional variable */
2421  sdpasol = sdpisolver->sdpa->getResultXVec();
2422  /* 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) */
2423  for (v = 0; v < sdpisolver->nvars; v++)
2424  {
2425  if ( sdpisolver->inputtosdpamapper[v] > 0 )
2426  {
2427  /* minus one because the inputtosdpamapper gives the sdpa indices which start at one, but the array starts at 0 */
2428  dualsol[v] = sdpasol[sdpisolver->inputtosdpamapper[v] - 1];
2429  }
2430  else
2431  {
2432  /* this is the value that was saved when inserting, as this variable has lb=ub */
2433  assert( -sdpisolver->inputtosdpamapper[v] <= sdpisolver->nvars - sdpisolver->nactivevars );
2434  dualsol[v] = sdpisolver->fixedvarsval[(-1 * sdpisolver->inputtosdpamapper[v]) - 1]; /*lint !e679*/
2435  }
2436  }
2437  }
2438  return SCIP_OKAY;
2439 }
2440 
2442 SCIP_RETCODE SCIPsdpiSolverGetPreoptimalPrimalNonzeros(
2443  SCIP_SDPISOLVER* sdpisolver,
2444  int nblocks,
2445  int* startXnblocknonz
2447  )
2448 {
2449  int b;
2450  int sdpablock;
2451  int sdpaind;
2452  int i;
2453  int c;
2454  int r;
2455  int blocksize;
2456 
2457  assert( sdpisolver != NULL );
2458  assert( nblocks > 0 );
2459  assert( startXnblocknonz != NULL );
2460  CHECK_IF_SOLVED( sdpisolver );
2461 
2462  if ( ! sdpisolver->preoptimalsolexists )
2463  {
2464  SCIPdebugMessage("Failed to retrieve preoptimal solution for warmstarting purposes. \n");
2465  startXnblocknonz[0] = -1;
2466  return SCIP_OKAY;
2467  }
2468 
2469  if ( nblocks != sdpisolver->nsdpblocks + 1 )
2470  {
2471  SCIPerrorMessage("SCIPsdpiSolverGetPreoptimalPrimalNonzeros expected nblocks = %d but got %d\n", sdpisolver->nsdpblocks + 1, nblocks);
2472  return SCIP_LPERROR;
2473  }
2474 
2475  /* iterate over all SDP-blocks, get the solution and count the nonzeros */
2476  for (b = 0; b < sdpisolver->nsdpblocks; b++)
2477  {
2478  sdpablock = sdpisolver->inputtoblockmapper[b];
2479 
2480  startXnblocknonz[b] = 0;
2481 
2482  if ( sdpablock != -1 )
2483  {
2484  blocksize = (int) sdpisolver->sdpa->getBlockSize(sdpablock);
2485 
2486  /* iterate once over the upper triangular part of the matrix (saving the corresponding entries in the lower triangular part for the SDPI) */
2487  for (r = 0; r < blocksize; r++)
2488  {
2489  for (c = r; c < blocksize; c++)
2490  {
2491  sdpaind = r + blocksize * c;
2492 
2493  if ( REALABS(sdpisolver->preoptimalsolx[sdpablock - 1][sdpaind]) > sdpisolver->epsilon )
2494  startXnblocknonz[b]++;
2495  }
2496  }
2497  }
2498  }
2499 
2500  /* compute the entry for the LP-block */
2501  startXnblocknonz[nblocks - 1] = 0;
2502  sdpablock = (int) sdpisolver->sdpa->getBlockNumber();
2503 
2504  if ( sdpisolver->sdpa->getBlockType(sdpablock) == SDPA::LP )
2505  {
2506  blocksize = (int) sdpisolver->sdpa->getBlockSize(sdpablock);
2507 
2508  for (i = 0; i < blocksize; i++)
2509  {
2510  if ( REALABS(sdpisolver->preoptimalsolxlp[i]) > sdpisolver->epsilon )
2511  startXnblocknonz[b]++;
2512  }
2513  }
2514 
2515  return SCIP_OKAY;
2516 }
2517 
2518 
2526 SCIP_RETCODE SCIPsdpiSolverGetPreoptimalSol(
2527  SCIP_SDPISOLVER* sdpisolver,
2528  SCIP_Bool* success,
2529  SCIP_Real* dualsol,
2530  int* dualsollength,
2532  int nblocks,
2533  int* startXnblocknonz,
2535  int** startXrow,
2536  int** startXcol,
2537  SCIP_Real** startXval
2538  )
2539 {
2540  int v;
2541  int b;
2542  int r;
2543  int c;
2544  int sdpaind;
2545  int sdpablock;
2546  int blocksize;
2547  int blocknnonz;
2548  SCIP_Bool msgthrown = FALSE;
2549 
2550  assert( sdpisolver != NULL );
2551  assert( success != NULL );
2552  assert( dualsol != NULL );
2553  assert( dualsollength != NULL );
2554  assert( *dualsollength >= 0 );
2555  assert( startXnblocknonz != NULL || nblocks == -1 );
2556  assert( startXrow != NULL || nblocks == -1 );
2557  assert( startXcol != NULL || nblocks == -1 );
2558  assert( startXval != NULL || nblocks == -1 );
2559 
2560  if ( ! sdpisolver->preoptimalsolexists )
2561  {
2562  SCIPdebugMessage("Failed to retrieve preoptimal solution for warmstarting purposes. \n");
2563  *success = FALSE;
2564  return SCIP_OKAY;
2565  }
2566 
2567  if ( *dualsollength < sdpisolver->nvars )
2568  {
2569  SCIPdebugMessage("Insufficient memory in SCIPsdpiSolverGetPreoptimalSol: needed %d, given %d\n", sdpisolver->nvars, *dualsollength);
2570  *success = FALSE;
2571  *dualsollength = sdpisolver->nvars;
2572  return SCIP_OKAY;
2573  }
2574 
2575  for (v = 0; v < sdpisolver->nvars; v++)
2576  {
2577  if (sdpisolver->inputtosdpamapper[v] > 0)
2578  {
2579  /* minus one because the inputtosdpamapper gives the dsdp indices which start at one, but the array starts at 0 */
2580  dualsol[v] = sdpisolver->preoptimalsol[sdpisolver->inputtosdpamapper[v] - 1];
2581  }
2582  else
2583  {
2584  /* this is the value that was saved when inserting, as this variable has lb=ub */
2585  dualsol[v] = sdpisolver->fixedvarsval[(-1 * sdpisolver->inputtosdpamapper[v]) - 1]; /*lint !e679*/
2586  }
2587  }
2588 
2589  *dualsollength = sdpisolver->nvars;
2590 
2591  if ( nblocks > -1 )
2592  {
2593  assert( startXnblocknonz != NULL );
2594  assert( startXrow != NULL );
2595  assert( startXcol != NULL );
2596  assert( startXval != NULL );
2597 
2598  /* iterate over all SDP-blocks and get the solution */
2599  for (b = 0; b < sdpisolver->nsdpblocks; b++)
2600  {
2601  assert( startXrow[b] != NULL );
2602  assert( startXcol[b] != NULL );
2603  assert( startXval[b] != NULL );
2604 
2605  sdpablock = sdpisolver->inputtoblockmapper[b];
2606 
2607  blocknnonz = 0;
2608 
2609  if ( sdpablock != -1 )
2610  {
2611  /* since we reset the preoptimalsolution for every solve, the blocksize should have stayed the same */
2612  blocksize = (int) sdpisolver->sdpa->getBlockSize(sdpablock);
2613  blocknnonz = 0;
2614 
2615  /* iterate once over the upper triangular part of the matrix (saving the corresponding entries in the lower triangular part for the SDPI) */
2616  for (r = 0; r < blocksize; r++)
2617  {
2618  for (c = r; c < blocksize; c++)
2619  {
2620  sdpaind = r + blocksize * c;
2621 
2622  if ( REALABS(sdpisolver->preoptimalsolx[sdpablock - 1][sdpaind]) > sdpisolver->epsilon )
2623  {
2624  if ( blocknnonz < startXnblocknonz[b] )
2625  {
2626  startXrow[b][blocknnonz] = sdpisolver->blockindmapper[b][c];
2627  startXcol[b][blocknnonz] = sdpisolver->blockindmapper[b][r];
2628  startXval[b][blocknnonz] = sdpisolver->preoptimalsolx[sdpablock - 1][sdpaind]; /* -1 because sdpa starts counting at 1 */
2629  blocknnonz++;
2630  }
2631  else
2632  {
2633  if ( ! msgthrown )
2634  {
2635  SCIPdebugMessage("Unsufficient arraylength %d for block %d in SCIPsdpiSolverGetPrimalMatrix!\n", startXnblocknonz[b], b);
2636  msgthrown = TRUE;
2637  }
2638  }
2639  }
2640  }
2641  }
2642 
2643  startXnblocknonz[b] = blocknnonz;
2644  }
2645  else
2646  startXnblocknonz[b] = 0;
2647  }
2648 
2649  /* compute entries for the LP-block */
2650  blocknnonz = 0;
2651 
2652  /* since we reset the preoptimalsolution for every solve, the number of blocks should have stayed the same */
2653  sdpablock = (int) sdpisolver->sdpa->getBlockNumber();
2654 
2655  if ( sdpisolver->sdpa->getBlockType(sdpablock) == SDPA::LP )
2656  {
2657  int i;
2658 
2659  /* since we reset the preoptimalsolution for every solve, the blocksize should have stayed the same */
2660  blocksize = (int) sdpisolver->sdpa->getBlockSize(sdpablock);
2661 
2662  /* iterate over LP constraints */
2663  for (i = 0; i < blocksize - sdpisolver->nvarbounds; i++)
2664  {
2665  if ( REALABS(sdpisolver->preoptimalsolxlp[i]) > sdpisolver->epsilon )
2666  {
2667  if ( blocknnonz < startXnblocknonz[b] )
2668  {
2669  startXrow[b][blocknnonz] = sdpisolver->rowtoinputmapper[i];
2670  startXcol[b][blocknnonz] = sdpisolver->rowtoinputmapper[i];
2671  startXval[b][blocknnonz] = sdpisolver->preoptimalsolxlp[i]; /* -1 because sdpa starts counting at 1 */
2672  blocknnonz++;
2673  }
2674  else
2675  {
2676  blocknnonz++;
2677  if ( ! msgthrown )
2678  {
2679  SCIPdebugMessage("Unsufficient arraylength %d for LP block in SCIPsdpiSolverGetPrimalMatrix!\n", startXnblocknonz[b]);
2680  msgthrown = TRUE;
2681  }
2682  }
2683  }
2684  }
2685 
2686  /* iterate over varbounds */
2687  for (i = blocksize - sdpisolver->nvarbounds; i < blocksize; i++)
2688  {
2689  if ( REALABS(sdpisolver->preoptimalsolxlp[i]) > sdpisolver->epsilon )
2690  {
2691  if ( blocknnonz < startXnblocknonz[b] )
2692  {
2693  int inputpos;
2694  int vbpos;
2695 
2696  vbpos = i - (blocksize - sdpisolver->nvarbounds); /* position in varboundpos array */
2697 
2698  /* inputpos is 2 * nlprows (for lhs and rhs) + 2 * inputvar for lb and 2 * inputvar + 1 for ub */
2699  if ( sdpisolver->varboundpos[vbpos] > 0 ) /* rhs */
2700  inputpos = 2 * sdpisolver->nsdpalpcons + 2 * sdpisolver->sdpatoinputmapper[sdpisolver->varboundpos[vbpos] - 1] + 1;
2701  else
2702  inputpos = 2 * sdpisolver->nsdpalpcons + 2 * sdpisolver->sdpatoinputmapper[-1 * sdpisolver->varboundpos[vbpos] - 1];
2703  startXrow[b][blocknnonz] = inputpos;
2704  startXcol[b][blocknnonz] = inputpos;
2705  startXval[b][blocknnonz] = sdpisolver->preoptimalsolxlp[i]; /* -1 because sdpa starts counting at 1 */
2706  blocknnonz++;
2707  }
2708  else
2709  {
2710  blocknnonz++;
2711  if ( ! msgthrown )
2712  {
2713  SCIPdebugMessage("Insufficient arraylength %d for LP block & varbounds in SCIPsdpiSolverGetPrimalMatrix!\n", startXnblocknonz[b]);
2714  msgthrown = TRUE;
2715  }
2716  }
2717  }
2718  }
2719  startXnblocknonz[b] = blocknnonz;
2720  }
2721  else
2722  startXnblocknonz[b] = 0;
2723  }
2724 
2725  *success = TRUE;
2726 
2727  return SCIP_OKAY;
2728 }
2729 
2737 SCIP_RETCODE SCIPsdpiSolverGetPrimalBoundVars(
2738  SCIP_SDPISOLVER* sdpisolver,
2739  SCIP_Real* lbvars,
2740  SCIP_Real* ubvars,
2741  int* arraylength
2743  )
2744 {/*lint !e1784*/
2745  int i;
2746  SCIP_Real* X; /* block of primal solution matrix corresponding to the LP-part */
2747  SDPA_INT lpblockind;
2748  int penoffset = 0;
2749  int nlpcons;
2750  int sdpapos;
2751  int pos;
2752 
2753  assert( sdpisolver != NULL );
2754  assert( sdpisolver->sdpa != NULL );
2755  assert( lbvars != NULL );
2756  assert( ubvars != NULL );
2757  assert( arraylength != NULL );
2758  assert( *arraylength >= 0 );
2759  CHECK_IF_SOLVED( sdpisolver );
2760 
2761  /* check if the arrays are long enough */
2762  if ( *arraylength < sdpisolver->nvars )
2763  {
2764  *arraylength = sdpisolver->nvars;
2765  SCIPdebugMessage("Insufficient length of array in SCIPsdpiSolverGetPrimalBoundVars (gave %d, needed %d)\n", *arraylength, sdpisolver->nvars);
2766  return SCIP_OKAY;
2767  }
2768 
2769  /* initialize the return-arrays with zero */
2770  for (i = 0; i < sdpisolver->nvars; i++)
2771  {
2772  lbvars[i] = 0.0;
2773  ubvars[i] = 0.0;
2774  }
2775 
2776  /* 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) */
2777  if ( sdpisolver->nvarbounds == 0 )
2778  {
2779  SCIPdebugMessage("Asked for PrimalBoundVars, but there were no variable bounds in sdpa, returning zero vector!\n");
2780  return SCIP_OKAY;
2781  }
2782 
2783  /* get the block of primal solution matrix corresponding to the LP-part from sdpa */
2784  lpblockind = (int) sdpisolver->sdpa->getBlockNumber(); /* the LP block is the last one and sdpa counts from one */
2785  if ( sdpisolver->sdpa->getBlockType(lpblockind) != SDPA::LP )
2786  {
2787  /* if the last block is not an LP-block, no variable bounds existed */
2788  *arraylength = 0;
2789 
2790  return SCIP_OKAY;
2791  }
2792  nlpcons = (int) sdpisolver->sdpa->getBlockSize(lpblockind);
2793  assert( nlpcons >= 0 );
2794 
2795  X = sdpisolver->sdpa->getResultYMat(lpblockind);
2796 
2797  /* iterate over all variable bounds and insert the corresponding primal variables in the right positions of the return-arrays */
2798  assert( sdpisolver->nvarbounds <= 2 * sdpisolver->nvars || (sdpisolver->nvarbounds <= 2 * sdpisolver->nvars + 1 && sdpisolver->penalty ) );
2799 
2800  /* if we solved a penalty formulation, the last variable bound belongs to the penalty variable, which we aren't interested in here */
2801  if ( sdpisolver->penalty )
2802  penoffset = 1;
2803  for (i = 0; i < sdpisolver->nvarbounds - penoffset; i++)
2804  {
2805  /* if it is a lower bound */
2806  if ( sdpisolver->varboundpos[i] < 0 )
2807  {
2808  sdpapos = -sdpisolver->varboundpos[i] - 1;
2809  assert( 0 <= sdpapos && sdpapos < sdpisolver->nactivevars );
2810  pos = sdpisolver->sdpatoinputmapper[sdpapos];
2811  assert( 0 <= pos && pos < sdpisolver->nvars );
2812  /* the last nvarbounds entries correspond to the varbounds */
2813  lbvars[pos] = X[nlpcons - sdpisolver->nvarbounds + i];
2814  }
2815  else /* if it is an upper bound */
2816  {
2817  sdpapos = sdpisolver->varboundpos[i] - 1;
2818  assert( 0 <= sdpapos && sdpapos < sdpisolver->nactivevars );
2819  pos = sdpisolver->sdpatoinputmapper[sdpapos];
2820  assert( 0 <= pos && pos < sdpisolver->nvars );
2821  /* the last nvarbounds entries correspond to the varbounds */
2822  ubvars[pos] = X[nlpcons - sdpisolver->nvarbounds + i]; /*lint !e679, !e834 */
2823  }
2824  }
2825 
2826  return SCIP_OKAY;
2827 }
2828 
2830 SCIP_RETCODE SCIPsdpiSolverGetPrimalNonzeros(
2831  SCIP_SDPISOLVER* sdpisolver,
2832  int nblocks,
2833  int* startXnblocknonz
2834  )
2835 {
2836  SCIP_Real* X;
2837  int b;
2838  int sdpablock;
2839  int sdpaind;
2840  int i;
2841  int r;
2842  int c;
2843  int blocksize;
2844 
2845  assert( sdpisolver != NULL );
2846  assert( nblocks > 0 );
2847  assert( startXnblocknonz != NULL );
2848  CHECK_IF_SOLVED( sdpisolver );
2849 
2850  if ( nblocks != sdpisolver->nsdpblocks + 1 )
2851  {
2852  SCIPerrorMessage("SCIPsdpiSolverGetPrimalNonzeros expected nblocks = %d but got %d\n", sdpisolver->nsdpblocks + 1, nblocks);
2853  return SCIP_LPERROR;
2854  }
2855 
2856  /* iterate over all SDP-blocks, get the solution and count the nonzeros */
2857  for (b = 0; b < sdpisolver->nsdpblocks; b++)
2858  {
2859  sdpablock = sdpisolver->inputtoblockmapper[b];
2860 
2861  startXnblocknonz[b] = 0;
2862 
2863  if ( sdpablock != -1 )
2864  {
2865  X = sdpisolver->sdpa->getResultXMat(sdpablock);
2866  blocksize = (int) sdpisolver->sdpa->getBlockSize(sdpablock);
2867 
2868  /* iterate once over the upper triangular part of the matrix (saving the corresponding entries in the lower triangular part for the SDPI) */
2869  for (r = 0; r < blocksize; r++)
2870  {
2871  for (c = r; c < blocksize; c++)
2872  {
2873  sdpaind = r + blocksize * c;
2874 
2875  if ( REALABS(X[sdpaind]) > sdpisolver->epsilon )
2876  startXnblocknonz[b]++;
2877  }
2878  }
2879  }
2880  }
2881 
2882  /* compute the entry for the LP-block */
2883  startXnblocknonz[nblocks - 1] = 0;
2884  sdpablock = (int) sdpisolver->sdpa->getBlockNumber();
2885 
2886  if ( sdpisolver->sdpa->getBlockType(sdpablock) == SDPA::LP )
2887  {
2888  X = sdpisolver->sdpa->getResultXMat(sdpablock);
2889  blocksize = (int) sdpisolver->sdpa->getBlockSize(sdpablock);
2890 
2891  for (i = 0; i < blocksize; i++)
2892  {
2893  if ( REALABS(X[i]) > sdpisolver->epsilon )
2894  startXnblocknonz[b]++;
2895  }
2896  }
2897 
2898  return SCIP_OKAY;
2899 }
2900 
2907 SCIP_RETCODE SCIPsdpiSolverGetPrimalMatrix(
2908  SCIP_SDPISOLVER* sdpisolver,
2909  int nblocks,
2910  int* startXnblocknonz,
2912  int** startXrow,
2913  int** startXcol,
2914  SCIP_Real** startXval
2915  )
2916 {
2917  SCIP_Real* X;
2918  int b;
2919  int r;
2920  int c;
2921  int sdpaind;
2922  int sdpablock;
2923  int blocksize;
2924  int blocknnonz;
2925  SCIP_Bool msgthrown = FALSE;
2926 
2927  assert( sdpisolver != NULL );
2928  assert( nblocks > 0 );
2929  assert( startXnblocknonz != NULL );
2930  assert( startXrow != NULL );
2931  assert( startXcol != NULL );
2932  assert( startXval != NULL );
2933  CHECK_IF_SOLVED( sdpisolver );
2934 
2935  if ( nblocks != sdpisolver->nsdpblocks + 1 )
2936  {
2937  SCIPerrorMessage("SCIPsdpiSolverGetPrimalNonzeros expected nblocks = %d but got %d\n", sdpisolver->nsdpblocks + 1, nblocks);
2938  return SCIP_LPERROR;
2939  }
2940 
2941  /* iterate over all SDP-blocks and get the solution */
2942  for (b = 0; b < sdpisolver->nsdpblocks; b++)
2943  {
2944  sdpablock = sdpisolver->inputtoblockmapper[b];
2945 
2946  blocknnonz = 0;
2947 
2948  if ( sdpablock != -1 )
2949  {
2950  X = sdpisolver->sdpa->getResultYMat(sdpablock);
2951  blocksize = (int) sdpisolver->sdpa->getBlockSize(sdpablock);
2952  blocknnonz = 0;
2953 
2954  /* iterate once over the upper triangular part of the matrix (saving the corresponding entries in the lower triangular part for the SDPI) */
2955  for (r = 0; r < blocksize; r++)
2956  {
2957  for (c = r; c < blocksize; c++)
2958  {
2959  sdpaind = r + blocksize * c;
2960 
2961  if ( REALABS(X[sdpaind]) > sdpisolver->epsilon )
2962  {
2963  if ( blocknnonz < startXnblocknonz[b] )
2964  {
2965  startXrow[b][blocknnonz] = sdpisolver->blockindmapper[b][c];
2966  startXcol[b][blocknnonz] = sdpisolver->blockindmapper[b][r];
2967  startXval[b][blocknnonz] = X[sdpaind];
2968  blocknnonz++;
2969  }
2970  else
2971  {
2972  if ( ! msgthrown )
2973  {
2974  SCIPdebugMessage("Insufficient arraylength %d for block %d in SCIPsdpiSolverGetPrimalMatrix!\n", startXnblocknonz[b], b);
2975  msgthrown = TRUE;
2976  }
2977  }
2978  }
2979  }
2980  }
2981 
2982  startXnblocknonz[b] = blocknnonz;
2983  }
2984  else
2985  startXnblocknonz[b] = 0;
2986  }
2987 
2988  /* compute entries for the LP-block */
2989  blocknnonz = 0;
2990  sdpablock = (int) sdpisolver->sdpa->getBlockNumber();
2991 
2992  if ( sdpisolver->sdpa->getBlockType(sdpablock) == SDPA::LP )
2993  {
2994  int i;
2995 
2996  X = sdpisolver->sdpa->getResultXMat(sdpablock);
2997  blocksize = (int) sdpisolver->sdpa->getBlockSize(sdpablock);
2998 
2999  /* iterate over LP constraints */
3000  for (i = 0; i < blocksize - sdpisolver->nvarbounds; i++)
3001  {
3002  if ( REALABS(X[i]) > sdpisolver->epsilon )
3003  {
3004  if ( blocknnonz < startXnblocknonz[b] )
3005  {
3006  startXrow[b][blocknnonz] = sdpisolver->rowtoinputmapper[i];
3007  startXcol[b][blocknnonz] = sdpisolver->rowtoinputmapper[i];
3008  startXval[b][blocknnonz] = X[i];
3009  blocknnonz++;
3010  }
3011  else
3012  {
3013  blocknnonz++;
3014  if ( ! msgthrown )
3015  {
3016  SCIPdebugMessage("Unsufficient arraylength %d for LP block in SCIPsdpiSolverGetPrimalMatrix!\n", startXnblocknonz[b]);
3017  msgthrown = TRUE;
3018  }
3019  }
3020  }
3021  }
3022 
3023  /* iterate over varbounds */
3024  for (i = blocksize - sdpisolver->nvarbounds; i < blocksize; i++)
3025  {
3026  if ( REALABS(X[i]) > sdpisolver->epsilon )
3027  {
3028  if ( blocknnonz < startXnblocknonz[b] )
3029  {
3030  int inputpos;
3031  int vbpos;
3032 
3033  vbpos = i - (blocksize - sdpisolver->nvarbounds); /* position in varboundpos array */
3034 
3035  /* inputpos is 2 * nlprows (for lhs and rhs) + 2 * inputvar for lb and 2 * inputvar + 1 for ub */
3036  if ( sdpisolver->varboundpos[vbpos] > 0 ) /* rhs */
3037  inputpos = 2 * sdpisolver->nsdpalpcons + 2 * sdpisolver->sdpatoinputmapper[sdpisolver->varboundpos[vbpos] - 1] + 1;
3038  else
3039  inputpos = 2 * sdpisolver->nsdpalpcons + 2 * sdpisolver->sdpatoinputmapper[-1 * sdpisolver->varboundpos[vbpos] - 1];
3040  startXrow[b][blocknnonz] = inputpos;
3041  startXcol[b][blocknnonz] = inputpos;
3042  startXval[b][blocknnonz] = X[i];
3043  blocknnonz++;
3044  }
3045  else
3046  {
3047  blocknnonz++;
3048  if ( ! msgthrown )
3049  {
3050  SCIPdebugMessage("Unsufficient arraylength %d for LP block & varbounds in SCIPsdpiSolverGetPrimalMatrix!\n", startXnblocknonz[b]);
3051  msgthrown = TRUE;
3052  }
3053  }
3054  }
3055  }
3056  startXnblocknonz[b] = blocknnonz;
3057  }
3058  else
3059  startXnblocknonz[b] = 0;
3060 
3061  return SCIP_OKAY;
3062 }
3063 
3065 SCIP_Real SCIPsdpiSolverGetMaxPrimalEntry(
3066  SCIP_SDPISOLVER* sdpisolver
3067  )
3068 {
3069  SCIP_Real* X;
3070  SCIP_Real maxentry;
3071  int nblocks;
3072  int blocksize;
3073  int arraylength;
3074  int b;
3075  int i;
3076 
3077  assert( sdpisolver != NULL );
3078  assert( sdpisolver->sdpa != NULL );
3079 
3080  nblocks = (int) sdpisolver->sdpa->getBlockNumber();
3081  maxentry = 0.0;
3082 
3083  /* iterate over all blocks */
3084  for (b = 0; b < nblocks; b++)
3085  {
3086  /* get the length of the X-array, if it is an LP block, the array has length n, otherwise n*(n+1)/2 */
3087  blocksize = (int) sdpisolver->sdpa->getBlockSize((long long) b + 1);
3088  arraylength = sdpisolver->sdpa->getBlockType((long long) b + 1) == SDPA::LP ? blocksize : (blocksize * (blocksize + 1) / 2);
3089 
3090  /* iterate over all entries in this block */
3091  X = sdpisolver->sdpa->getResultYMat((long long) b + 1);/*lint !e747*/
3092  for (i = 0; i < arraylength; i++)
3093  {
3094  if ( X[i] > maxentry )
3095  maxentry = X[i];
3096  }
3097  }
3098 
3099  return maxentry;
3100 }
3101 
3103 SCIP_RETCODE SCIPsdpiSolverGetTime(
3104  SCIP_SDPISOLVER* sdpisolver,
3105  SCIP_Real* opttime
3106  )
3107 {
3108  assert( sdpisolver != NULL );
3109  assert( opttime != NULL );
3110 
3111  /* not implemented */
3112  *opttime = 0.0;
3113 
3114  return SCIP_OKAY;
3115 }
3116 
3118 SCIP_RETCODE SCIPsdpiSolverGetIterations(
3119  SCIP_SDPISOLVER* sdpisolver,
3120  int* iterations
3121  )
3122 {/*lint !e1784*/
3123  assert( sdpisolver != NULL );
3124  assert( sdpisolver->sdpa != NULL );
3125  assert( iterations != NULL );
3126 
3127  *iterations = sdpisolver->niterations;
3128 
3129  return SCIP_OKAY;
3130 }
3131 
3133 SCIP_RETCODE SCIPsdpiSolverGetSdpCalls(
3134  SCIP_SDPISOLVER* sdpisolver,
3135  int* calls
3136  )
3137 {/*lint !e1784*/
3138  assert( sdpisolver != NULL );
3139  assert( sdpisolver->sdpa != NULL );
3140  assert( calls != NULL );
3141 
3142  *calls = sdpisolver->nsdpcalls;
3143 
3144  return SCIP_OKAY;
3145 }
3146 
3148 SCIP_RETCODE SCIPsdpiSolverSettingsUsed(
3149  SCIP_SDPISOLVER* sdpisolver,
3150  SCIP_SDPSOLVERSETTING* usedsetting
3151  )
3152 {/*lint !e1784*/
3153  assert( sdpisolver != NULL );
3154  assert( usedsetting != NULL );
3155  CHECK_IF_SOLVED(sdpisolver);
3156 
3157  *usedsetting = sdpisolver->usedsetting;
3158 
3159  return SCIP_OKAY;
3160 }
3161 
3167 /*
3168  * Numerical Methods
3169  */
3170 
3175 SCIP_Real SCIPsdpiSolverInfinity(
3176  SCIP_SDPISOLVER* sdpisolver
3177  )
3178 { /*lint --e{715}*/
3179  return 1E+20; /* default infinity from SCIP */
3180 }
3181 
3183 SCIP_Bool SCIPsdpiSolverIsInfinity(
3184  SCIP_SDPISOLVER* sdpisolver,
3185  SCIP_Real val
3186  )
3187 {/*lint !e1784*/
3188  return ( val <= -SCIPsdpiSolverInfinity(sdpisolver) || val >= SCIPsdpiSolverInfinity(sdpisolver) );
3189 }
3190 
3192 SCIP_RETCODE SCIPsdpiSolverGetRealpar(
3193  SCIP_SDPISOLVER* sdpisolver,
3194  SCIP_SDPPARAM type,
3195  SCIP_Real* dval
3196  )
3197 {
3198  assert( sdpisolver != NULL );
3199  assert( dval != NULL );
3200 
3201  switch( type )/*lint --e{788}*/
3202  {
3203  case SCIP_SDPPAR_EPSILON:
3204  *dval = sdpisolver->epsilon;
3205  break;
3206  case SCIP_SDPPAR_GAPTOL:
3207  *dval = sdpisolver->gaptol;
3208  break;
3209  case SCIP_SDPPAR_FEASTOL:
3210  *dval = sdpisolver->feastol;
3211  break;
3212  case SCIP_SDPPAR_SDPSOLVERFEASTOL:
3213  *dval = sdpisolver->sdpsolverfeastol;
3214  break;
3215  case SCIP_SDPPAR_PENALTYPARAM:
3216  *dval = 0.0;
3217  SCIPdebugMessage("Parameter SCIP_SDPPAR_PENALTYPARAM not used by SDPA"); /* this parameter is only used by DSDP */
3218  break;
3219  case SCIP_SDPPAR_OBJLIMIT:
3220  *dval = sdpisolver->objlimit;
3221  break;
3222  case SCIP_SDPPAR_LAMBDASTAR:
3223  *dval = sdpisolver->lambdastar;
3224  break;
3225  case SCIP_SDPPAR_WARMSTARTPOGAP:
3226  *dval = sdpisolver->preoptimalgap;
3227  break;
3228  default:
3229  return SCIP_PARAMETERUNKNOWN;
3230  }
3231 
3232  return SCIP_OKAY;
3233 }
3234 
3236 SCIP_RETCODE SCIPsdpiSolverSetRealpar(
3237  SCIP_SDPISOLVER* sdpisolver,
3238  SCIP_SDPPARAM type,
3239  SCIP_Real dval
3240  )
3241 {/*lint !e1784*/
3242  assert( sdpisolver != NULL );
3243 
3244  switch( type )/*lint --e{788}*/
3245  {
3246  case SCIP_SDPPAR_EPSILON:
3247  sdpisolver->epsilon = dval;
3248  SCIPdebugMessage("Setting sdpisolver epsilon to %g.\n", dval);
3249  break;
3250  case SCIP_SDPPAR_GAPTOL:
3251  sdpisolver->gaptol = dval;
3252  SCIPdebugMessage("Setting sdpisolver gaptol to %g.\n", dval);
3253  break;
3254  case SCIP_SDPPAR_FEASTOL:
3255  sdpisolver->feastol = dval;
3256  SCIPdebugMessage("Setting sdpisolver feastol to %g.\n", dval);
3257  break;
3258  case SCIP_SDPPAR_SDPSOLVERFEASTOL:
3259  sdpisolver->sdpsolverfeastol = dval;
3260  SCIPdebugMessage("Setting sdpisolver sdpsolverfeastol to %g.\n", dval);
3261  break;
3262  case SCIP_SDPPAR_PENALTYPARAM:
3263  SCIPdebugMessage("Parameter SCIP_SDPPAR_PENALTYPARAM not used by SDPA\n."); /* this parameter is only used by DSDP */
3264  break;
3265  case SCIP_SDPPAR_OBJLIMIT:
3266  SCIPdebugMessage("Setting sdpisolver objlimit to %g.\n", dval);
3267  sdpisolver->objlimit = dval;
3268  break;
3269  case SCIP_SDPPAR_LAMBDASTAR:
3270  SCIPdebugMessage("Setting sdpisolver lambdastar parameter to %g.\n", dval);
3271  sdpisolver->lambdastar = dval;
3272  break;
3273  case SCIP_SDPPAR_WARMSTARTPOGAP:
3274  SCIPdebugMessage("Setting sdpisolver preoptgap to %g.\n", dval);
3275  sdpisolver->preoptimalgap = dval;
3276  break;
3277  default:
3278  return SCIP_PARAMETERUNKNOWN;
3279  }
3280 
3281  return SCIP_OKAY;
3282 }
3283 
3285 SCIP_RETCODE SCIPsdpiSolverGetIntpar(
3286  SCIP_SDPISOLVER* sdpisolver,
3287  SCIP_SDPPARAM type,
3288  int* ival
3289  )
3290 {/*lint !e1784*/
3291  assert( sdpisolver != NULL );
3292 
3293  switch( type )/*lint --e{788}*/
3294  {
3295  case SCIP_SDPPAR_SDPINFO:
3296  *ival = (int) sdpisolver->sdpinfo;
3297  SCIPdebugMessage("Getting sdpisolver information output (%d).\n", *ival);
3298  break;
3299  case SCIP_SDPPAR_NTHREADS:
3300  *ival = sdpisolver->nthreads;
3301  SCIPdebugMessage("Getting sdpisolver number of threads: %d.\n", *ival);
3302  break;
3303  default:
3304  return SCIP_PARAMETERUNKNOWN;
3305  }
3306 
3307  return SCIP_OKAY;
3308 }
3309 
3311 SCIP_RETCODE SCIPsdpiSolverSetIntpar(
3312  SCIP_SDPISOLVER* sdpisolver,
3313  SCIP_SDPPARAM type,
3314  int ival
3315  )
3316 {/*lint !e1784*/
3317  assert( sdpisolver != NULL );
3318 
3319  switch( type )/*lint --e{788}*/
3320  {
3321  case SCIP_SDPPAR_SDPINFO:
3322  sdpisolver->sdpinfo = (SCIP_Bool) ival;
3323  SCIPdebugMessage("Setting sdpisolver information output (%d).\n", ival);
3324  break;
3325  case SCIP_SDPPAR_NTHREADS:
3326  sdpisolver->nthreads = ival;
3327  SCIPdebugMessage("Setting sdpisolver number of threads to %d.\n", ival);
3328  break;
3329  default:
3330  return SCIP_PARAMETERUNKNOWN;
3331  }
3332 
3333  return SCIP_OKAY;
3334 }
3335 
3337 SCIP_RETCODE SCIPsdpiSolverComputeLambdastar(
3338  SCIP_SDPISOLVER* sdpisolver,
3339  SCIP_Real maxguess
3340  )
3341 {/*lint !e1784*/
3342  SCIP_Real compval;
3343 
3344  assert( sdpisolver != NULL );
3345 
3346  /* we set the value to min{max{MIN_LAMBDASTAR, LAMBDASTAR_FACTOR * MAX_GUESS}, MAX_LAMBDASTAR}, where MAX_GUESS is the maximum of the guesses
3347  * of the SDP-Blocks, if the define LAMBDASTAR_TWOPOINTS is set, we instead choose either LAMBDASTAR_LOW or HIGH depending on LAMBASTAR_THRESHOLD */
3348 
3349 #ifdef LAMBDASTAR_TWOPOINTS
3350  if ( maxguess < LAMBDASTAR_THRESHOLD )
3351  compval = LAMBDASTAR_LOW;
3352  else
3353  compval = LAMBDASTAR_HIGH;
3354 #else
3355  compval = LAMBDASTAR_FACTOR * maxguess;
3356 #endif
3357 
3358  if ( compval < MIN_LAMBDASTAR )
3359  {
3360  sdpisolver->lambdastar = MIN_LAMBDASTAR;
3361  SCIPdebugMessage("Setting lambdastar to %g.\n", MIN_LAMBDASTAR);
3362  }
3363  else if ( compval > MAX_LAMBDASTAR )
3364  {
3365  sdpisolver->lambdastar = MAX_LAMBDASTAR;
3366  SCIPdebugMessage("Setting lambdastar to %g.\n", MAX_LAMBDASTAR);
3367  }
3368  else
3369  {
3370  sdpisolver->lambdastar = compval;
3371  SCIPdebugMessage("Setting lambdastar to %g.\n", compval);
3372  }
3373 
3374  return SCIP_OKAY;
3375 }
3376 
3378 SCIP_RETCODE SCIPsdpiSolverComputePenaltyparam(
3379  SCIP_SDPISOLVER* sdpisolver,
3380  SCIP_Real maxcoeff,
3381  SCIP_Real* penaltyparam
3382  )
3383 {/*lint !e1784*/
3384  SCIP_Real compval;
3385 
3386  assert( sdpisolver != NULL );
3387  assert( penaltyparam != NULL );
3388 
3389  compval = PENALTYPARAM_FACTOR * maxcoeff;
3390 
3391  if ( compval < MIN_PENALTYPARAM )
3392  {
3393  SCIPdebugMessage("Setting penalty parameter to %g.\n", MIN_PENALTYPARAM);
3394  *penaltyparam = MIN_PENALTYPARAM;
3395  }
3396  else if ( compval > MAX_PENALTYPARAM )
3397  {
3398  SCIPdebugMessage("Setting penalty parameter to %g.\n", MAX_PENALTYPARAM);
3399  *penaltyparam = MAX_PENALTYPARAM;
3400  }
3401  else
3402  {
3403  SCIPdebugMessage("Setting penalty parameter to %g.\n", compval);
3404  *penaltyparam = compval;
3405  }
3406  return SCIP_OKAY;
3407 }
3408 
3410 SCIP_RETCODE SCIPsdpiSolverComputeMaxPenaltyparam(
3411  SCIP_SDPISOLVER* sdpisolver,
3412  SCIP_Real penaltyparam,
3413  SCIP_Real* maxpenaltyparam
3414  )
3415 {/*lint !e1784*/
3416  SCIP_Real compval;
3417 
3418  assert( sdpisolver != NULL );
3419  assert( maxpenaltyparam != NULL );
3420 
3421  compval = penaltyparam * MAXPENALTYPARAM_FACTOR;
3422 
3423  if ( compval < MAX_MAXPENALTYPARAM )
3424  {
3425  *maxpenaltyparam = compval;
3426  SCIPdebugMessage("Setting maximum penalty parameter to %g.\n", compval);
3427  }
3428  else
3429  {
3430  *maxpenaltyparam = MAX_MAXPENALTYPARAM;
3431  SCIPdebugMessage("Setting penalty parameter to %g.\n", MAX_MAXPENALTYPARAM);
3432  }
3433 
3434  return SCIP_OKAY;
3435 }
3436 
3442 /*
3443  * File Interface Methods
3444  */
3445 
3450 SCIP_RETCODE SCIPsdpiSolverReadSDP(
3451  SCIP_SDPISOLVER* sdpisolver,
3452  const char* fname
3453  )
3454 { /*lint --e{715}*/
3455  SCIPdebugMessage("Not implemented yet\n");
3456  return SCIP_LPERROR;
3457 }
3458 
3460 SCIP_RETCODE SCIPsdpiSolverWriteSDP(
3461  SCIP_SDPISOLVER* sdpisolver,
3462  const char* fname
3463  )
3464 {/*lint !e1784*/
3465  assert( fname != NULL );
3466 
3467  sdpisolver->sdpa->writeInputSparse(const_cast<char*>(fname), const_cast<char*>("%8.3f"));
3468 
3469  return SCIP_OKAY;
3470 }
3471 
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, SCIP_Real *lplhs, SCIP_Real *lprhs, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *starty, int *startZnblocknonz, int **startZrow, int **startZcol, SCIP_Real **startZval, int *startXnblocknonz, int **startXrow, int **startXcol, SCIP_Real **startXval, SCIP_SDPSOLVERSETTING startsettings, SCIP_Real timelimit, SDPI_CLOCK *usedsdpitime)
Definition: sdpisolver_sdpa.cpp:721
SCIP_SDPSOLVERSETTING_MEDIUM
Definition: type_sdpi.h:80
SCIPsdpiSolverGetSolverName
const char * SCIPsdpiSolverGetSolverName(void)
Definition: sdpisolver_sdpa.cpp:496
MIN_PENALTYPARAM
#define MIN_PENALTYPARAM
Definition: sdpisolver_sdpa.cpp:89
FEASTOLCHANGE
#define FEASTOLCHANGE
Definition: sdpisolver_sdpa.cpp:72
SCIP_SDPSOLVERSETTING_UNSOLVED
Definition: type_sdpi.h:77
SCIPsdpiSolverSetRealpar
SCIP_RETCODE SCIPsdpiSolverSetRealpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, SCIP_Real dval)
Definition: sdpisolver_sdpa.cpp:3237
SCIPsdpiSolverFree
SCIP_RETCODE SCIPsdpiSolverFree(SCIP_SDPISOLVER **sdpisolver)
Definition: sdpisolver_sdpa.cpp:621
SCIPsdpiSolverInfinity
SCIP_Real SCIPsdpiSolverInfinity(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:3176
SCIPsdpiSolverGetSolverPointer
void * SCIPsdpiSolverGetSolverPointer(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:517
SCIP_SDPPAR_EPSILON
Definition: type_sdpi.h:53
LAMBDASTAR_LOW
#define LAMBDASTAR_LOW
Definition: sdpisolver_sdpa.cpp:85
LAMBDASTAR_HIGH
#define LAMBDASTAR_HIGH
Definition: sdpisolver_sdpa.cpp:86
INFEASFEASTOLCHANGE
#define INFEASFEASTOLCHANGE
Definition: sdpisolver_sdpa.cpp:77
ensureMappingDataMemory
static SCIP_RETCODE ensureMappingDataMemory(SCIP_SDPISOLVER *sdpisolver, int nvars, int nsdpblocks, int *sdpblocksizes, int nlpcons, SCIP_Bool usepreoptimalsol)
Definition: sdpisolver_sdpa.cpp:247
SCIPsdpiSolverIgnoreInstability
SCIP_RETCODE SCIPsdpiSolverIgnoreInstability(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *success)
Definition: sdpisolver_sdpa.cpp:2300
SCIP_SDPSOLVERSETTING
enum SCIP_SDPSolverSetting SCIP_SDPSOLVERSETTING
Definition: type_sdpi.h:83
SCIPsdpiSolverIsPrimalUnbounded
SCIP_Bool SCIPsdpiSolverIsPrimalUnbounded(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1960
SCIPsdpSolcheckerCheck
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, SCIP_Real *lplhs, SCIP_Real *lprhs, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *solvector, SCIP_Real feastol, SCIP_Real epsilon, SCIP_Bool *infeasible)
Definition: sdpsolchecker.c:62
MAX_LAMBDASTAR
#define MAX_LAMBDASTAR
Definition: sdpisolver_sdpa.cpp:81
SCIPsdpiSolverGetObjval
SCIP_RETCODE SCIPsdpiSolverGetObjval(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *objval)
Definition: sdpisolver_sdpa.cpp:2311
PENALTYPARAM_FACTOR
#define PENALTYPARAM_FACTOR
Definition: sdpisolver_sdpa.cpp:91
SCIPsdpiSolverGetSolverDesc
const char * SCIPsdpiSolverGetSolverDesc(void)
Definition: sdpisolver_sdpa.cpp:504
SCIP_SDPSOLVERSETTING_FAST
Definition: type_sdpi.h:79
SCIPsdpiSolverIsInfinity
SCIP_Bool SCIPsdpiSolverIsInfinity(SCIP_SDPISOLVER *sdpisolver, SCIP_Real val)
Definition: sdpisolver_sdpa.cpp:3184
sdpisolver.h
interface methods for specific SDP-solvers
SCIP_SDPSOLVERSETTING_PENALTY
Definition: type_sdpi.h:78
SCIPsdpiSolverResetCounter
SCIP_RETCODE SCIPsdpiSolverResetCounter(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:684
SCIPsdpiSolverGetRealpar
SCIP_RETCODE SCIPsdpiSolverGetRealpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, SCIP_Real *dval)
Definition: sdpisolver_sdpa.cpp:3193
SCIPsdpiSolverIsObjlimExc
SCIP_Bool SCIPsdpiSolverIsObjlimExc(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2163
LAMBDASTAR_DEFAULT
#define LAMBDASTAR_DEFAULT
Definition: sdpisolver_sdpa.cpp:87
SCIPsdpiSolverGetSol
SCIP_RETCODE SCIPsdpiSolverGetSol(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *objval, SCIP_Real *dualsol, int *dualsollength)
Definition: sdpisolver_sdpa.cpp:2360
SDPI_Clock
Definition: struct_sdpiclock.h:72
MAXPENALTYPARAM_FACTOR
#define MAXPENALTYPARAM_FACTOR
Definition: sdpisolver_sdpa.cpp:93
SCIP_SDPPAR_PENALTYPARAM
Definition: type_sdpi.h:61
SCIPsdpiSolverGetPrimalMatrix
SCIP_RETCODE SCIPsdpiSolverGetPrimalMatrix(SCIP_SDPISOLVER *sdpisolver, int nblocks, int *startXnblocknonz, int **startXrow, int **startXcol, SCIP_Real **startXval)
Definition: sdpisolver_sdpa.cpp:2908
SCIP_SDPPAR_GAPTOL
Definition: type_sdpi.h:54
SCIPsdpiSolverFeasibilityKnown
SCIP_Bool SCIPsdpiSolverFeasibilityKnown(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1886
calcGrowSize
static int calcGrowSize(int initsize, int num)
Definition: sdpisolver_sdpa.cpp:212
MAX_MAXPENALTYPARAM
#define MAX_MAXPENALTYPARAM
Definition: sdpisolver_sdpa.cpp:92
SCIP_SDPPAR_OBJLIMIT
Definition: type_sdpi.h:57
SCIPsdpiSolverComputeLambdastar
SCIP_RETCODE SCIPsdpiSolverComputeLambdastar(SCIP_SDPISOLVER *sdpisolver, SCIP_Real maxguess)
Definition: sdpisolver_sdpa.cpp:3338
SCIPsdpiSolverIsPrimalFeasible
SCIP_Bool SCIPsdpiSolverIsPrimalFeasible(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2022
LAMBDASTAR_THRESHOLD
#define LAMBDASTAR_THRESHOLD
Definition: sdpisolver_sdpa.cpp:84
SCIPsdpiSolverDoesWarmstartNeedPrimal
SCIP_Bool SCIPsdpiSolverDoesWarmstartNeedPrimal(void)
Definition: sdpisolver_sdpa.cpp:534
SCIPsdpiSolverSettingsUsed
SCIP_RETCODE SCIPsdpiSolverSettingsUsed(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPSOLVERSETTING *usedsetting)
Definition: sdpisolver_sdpa.cpp:3149
SCIP_SDPPAR_WARMSTARTPOGAP
Definition: type_sdpi.h:66
SCIPsdpiSolverGetIntpar
SCIP_RETCODE SCIPsdpiSolverGetIntpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, int *ival)
Definition: sdpisolver_sdpa.cpp:3286
SCIP_SDPPAR_NTHREADS
Definition: type_sdpi.h:65
SCIPsdpiSolverWasSolved
SCIP_Bool SCIPsdpiSolverWasSolved(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1872
SCIPsdpiSolverCreate
SCIP_RETCODE SCIPsdpiSolverCreate(SCIP_SDPISOLVER **sdpisolver, SCIP_MESSAGEHDLR *messagehdlr, BMS_BLKMEM *blkmem, BMS_BUFMEM *bufmem)
Definition: sdpisolver_sdpa.cpp:552
LAMBDASTAR_FACTOR
#define LAMBDASTAR_FACTOR
Definition: sdpisolver_sdpa.cpp:82
CHECK_IF_SOLVED_BOOL
#define CHECK_IF_SOLVED_BOOL(sdpisolver)
Definition: sdpisolver_sdpa.cpp:119
sdpsolchecker.h
checks a given SDP solution for feasibility
checkFeastolAndResolve
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, SCIP_Real *lplhs, SCIP_Real *lprhs, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *feastol, SCIP_Real *gaptol)
Definition: sdpisolver_sdpa.cpp:354
isFixed
static SCIP_Bool isFixed(SCIP_SDPISOLVER *sdpisolver, SCIP_Real lb, SCIP_Real ub)
Definition: sdpisolver_sdpa.cpp:194
SCIP_SDPPAR_SDPINFO
Definition: type_sdpi.h:59
SCIPsdpiSolverGetPreoptimalSol
SCIP_RETCODE SCIPsdpiSolverGetPreoptimalSol(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *success, SCIP_Real *dualsol, int *dualsollength, int nblocks, int *startXnblocknonz, int **startXrow, int **startXcol, SCIP_Real **startXval)
Definition: sdpisolver_sdpa.cpp:2527
GAPTOLCHANGE
#define GAPTOLCHANGE
Definition: sdpisolver_sdpa.cpp:71
SCIPsdpiSolverGetDefaultSdpiSolverNpenaltyIncreases
int SCIPsdpiSolverGetDefaultSdpiSolverNpenaltyIncreases(void)
Definition: sdpisolver_sdpa.cpp:526
MIN_LAMBDASTAR
#define MIN_LAMBDASTAR
Definition: sdpisolver_sdpa.cpp:80
SCIP_SDPPAR_LAMBDASTAR
Definition: type_sdpi.h:64
SCIPsdpiSolverIncreaseCounter
SCIP_RETCODE SCIPsdpiSolverIncreaseCounter(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:672
SCIPsdpiSolverGetSolFeasibility
SCIP_RETCODE SCIPsdpiSolverGetSolFeasibility(SCIP_SDPISOLVER *sdpisolver, SCIP_Bool *primalfeasible, SCIP_Bool *dualfeasible)
Definition: sdpisolver_sdpa.cpp:1905
PENALTYBOUNDTOL
#define PENALTYBOUNDTOL
Definition: sdpisolver_sdpa.cpp:73
SCIPsdpiSolverWriteSDP
SCIP_RETCODE SCIPsdpiSolverWriteSDP(SCIP_SDPISOLVER *sdpisolver, const char *fname)
Definition: sdpisolver_sdpa.cpp:3461
SCIP_SDPSOLVERSETTING_STABLE
Definition: type_sdpi.h:81
SCIP_SDPPAR_FEASTOL
Definition: type_sdpi.h:55
SCIPsdpiSolverIsAcceptable
SCIP_Bool SCIPsdpiSolverIsAcceptable(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2274
SCIPsdpiSolverIsOptimal
SCIP_Bool SCIPsdpiSolverIsOptimal(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2251
CHECK_IF_SOLVED
#define CHECK_IF_SOLVED(sdpisolver)
Definition: sdpisolver_sdpa.cpp:108
SCIPsdpiSolverGetSdpCalls
SCIP_RETCODE SCIPsdpiSolverGetSdpCalls(SCIP_SDPISOLVER *sdpisolver, int *calls)
Definition: sdpisolver_sdpa.cpp:3134
SCIPsdpiSolverComputeMaxPenaltyparam
SCIP_RETCODE SCIPsdpiSolverComputeMaxPenaltyparam(SCIP_SDPISOLVER *sdpisolver, SCIP_Real penaltyparam, SCIP_Real *maxpenaltyparam)
Definition: sdpisolver_sdpa.cpp:3411
SCIPsdpiSolverGetTime
SCIP_RETCODE SCIPsdpiSolverGetTime(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *opttime)
Definition: sdpisolver_sdpa.cpp:3104
SCIPsdpiSolverGetPrimalBoundVars
SCIP_RETCODE SCIPsdpiSolverGetPrimalBoundVars(SCIP_SDPISOLVER *sdpisolver, SCIP_Real *lbvars, SCIP_Real *ubvars, int *arraylength)
Definition: sdpisolver_sdpa.cpp:2738
SDPIclockGetTime
SCIP_Real SDPIclockGetTime(SDPI_CLOCK *clck)
Definition: sdpiclock.c:250
SCIPsdpiSolverIsDualInfeasible
SCIP_Bool SCIPsdpiSolverIsDualInfeasible(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2084
SCIP_SDPPAR_SDPSOLVERFEASTOL
Definition: type_sdpi.h:56
SCIPsdpiSolverSetIntpar
SCIP_RETCODE SCIPsdpiSolverSetIntpar(SCIP_SDPISOLVER *sdpisolver, SCIP_SDPPARAM type, int ival)
Definition: sdpisolver_sdpa.cpp:3312
SCIPsdpiSolverIsPrimalInfeasible
SCIP_Bool SCIPsdpiSolverIsPrimalInfeasible(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:1991
SCIPsdpiSolverGetPrimalNonzeros
SCIP_RETCODE SCIPsdpiSolverGetPrimalNonzeros(SCIP_SDPISOLVER *sdpisolver, int nblocks, int *startXnblocknonz)
Definition: sdpisolver_sdpa.cpp:2831
SCIPsdpiSolverReadSDP
SCIP_RETCODE SCIPsdpiSolverReadSDP(SCIP_SDPISOLVER *sdpisolver, const char *fname)
Definition: sdpisolver_sdpa.cpp:3451
SCIP_SDPISOLVER
struct SCIP_SDPiSolver SCIP_SDPISOLVER
Definition: sdpisolver.h:71
SCIPsdpiSolverComputePenaltyparam
SCIP_RETCODE SCIPsdpiSolverComputePenaltyparam(SCIP_SDPISOLVER *sdpisolver, SCIP_Real maxcoeff, SCIP_Real *penaltyparam)
Definition: sdpisolver_sdpa.cpp:3379
SCIPsdpiSolverIsIterlimExc
SCIP_Bool SCIPsdpiSolverIsIterlimExc(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2182
SCIPsdpiSolverGetInternalStatus
int SCIPsdpiSolverGetInternalStatus(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2224
SCIP_SDPPARAM
enum SCIP_SDPParam SCIP_SDPPARAM
Definition: type_sdpi.h:72
SCIPsdpiSolverIsDualUnbounded
SCIP_Bool SCIPsdpiSolverIsDualUnbounded(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2053
INFEASMINFEASTOL
#define INFEASMINFEASTOL
Definition: sdpisolver_sdpa.cpp:78
SCIPsdpiSolverIsTimelimExc
SCIP_Bool SCIPsdpiSolverIsTimelimExc(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2204
SCIPsdpiSolverGetMaxPrimalEntry
SCIP_Real SCIPsdpiSolverGetMaxPrimalEntry(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:3066
SCIPsdpiSolverGetPreoptimalPrimalNonzeros
SCIP_RETCODE SCIPsdpiSolverGetPreoptimalPrimalNonzeros(SCIP_SDPISOLVER *sdpisolver, int nblocks, int *startXnblocknonz)
Definition: sdpisolver_sdpa.cpp:2443
BMS_CALL
#define BMS_CALL(x)
Definition: sdpisolver_sdpa.cpp:97
SCIPsdpiSolverIsDualFeasible
SCIP_Bool SCIPsdpiSolverIsDualFeasible(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2115
SCIPsdpiSolverGetIterations
SCIP_RETCODE SCIPsdpiSolverGetIterations(SCIP_SDPISOLVER *sdpisolver, int *iterations)
Definition: sdpisolver_sdpa.cpp:3119
SCIPsdpiSolverIsConverged
SCIP_Bool SCIPsdpiSolverIsConverged(SCIP_SDPISOLVER *sdpisolver)
Definition: sdpisolver_sdpa.cpp:2144
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, SCIP_Real *lplhs, SCIP_Real *lprhs, int lpnnonz, int *lprow, int *lpcol, SCIP_Real *lpval, SCIP_Real *starty, int *startZnblocknonz, int **startZrow, int **startZcol, SCIP_Real **startZval, int *startXnblocknonz, int **startXrow, int **startXcol, SCIP_Real **startXval, SCIP_SDPSOLVERSETTING startsettings, SCIP_Real timelimit, SDPI_CLOCK *usedsdpitime, SCIP_Bool *feasorig, SCIP_Bool *penaltybound)
Definition: sdpisolver_sdpa.cpp:808
MAX_PENALTYPARAM
#define MAX_PENALTYPARAM
Definition: sdpisolver_sdpa.cpp:90
Generated by   doxygen 1.8.11