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